A MILTISCALK SPLINE DERIVATIVE-BASF!) TRANSFORM
FOR IMAGE FUSION AND ENHANCEMENT
By
IZTOK KOREN
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1996
Copyright 1996
by
Iztok Koren
Mojim starsem
ACKNOWLEDGEMENTS
I would like to thank my mentors Dr. Fred J. Taylor and Dr. Andrew F. Laine
for their help and support. They gave me the freedom to explore research directions
that interested me most and provided me with the excellent work environment.
I am grateful to Dr. Jose C. Principe, Dr. John M. M. Anderson, and Dr.
Kermit Sigmon for serving on my committee and for their comments on my disser-
tation.
I would also like to thank all the members of the High Speed Digital Architec-
ture Laboratory and the Image Processing Research Group for their friendship and
help.
Last but not least, I want to thank my family for their support and under-
standing, and for bearing with me when I was unbearable.
IV
TABLE OF CONTENTS
ACKNOWLEDGEMENTS iv
LIST OF TABLES vii
LIST OF FIGURES viii
ABSTRACT ix
CHAPTERS
1 INTRODUCTION 1
1 . 1 Shortcomings of Traditional Methods of Wavelet Analysis .... 1
1.2 Thesis Overview 4
1.3 Notation 5
2 SIGNAL PROCESSING USING CENTRAL B-SPLINES 8
2.1 Central B-Splines: Definition and Properties 8
2.2 B-Spline Signal Interpolations 9
2.3 B-Spline Signal Approximations 14
3 TRANSFORM IN ONE DIMENSION 17
3.1 1-D Discrete Dyadic Wavelet Transform Revisited 17
3.2 Remarks 28
4 TRANSFORMS IN TWO DIMENSIONS 30
4.1 2-D Discrete Dyadic Wavelet Transform Revisited 30
4.2 Steerable Functions 36
4.3 Steerable Dyadic Wavelet Transform 38
4.4 Multiscale Spline Derivative-Based Transform 40
5 IMPLEMENTATION ISSUES 49
5.1 Finite Impulse Response Filters 49
5.2 Infinite Impulse Response Filters 55
6 APPLICATIONS 59
6.1 Image Fusion 59
6.2 Image Enhancement 61
7 CONCLUSION 66
REFERENCES 68
BIOGRAPHICAL SKETCH 72
VI
LIST OF TABLES
2.1 Transfer functions of direct B-spline filters for orders from to 9. . . 12
3.1 Impulse responses h(n), g(n), l(n), and k(n) for p = and d € {1,2,3}. 25
3.2 Impulse responses h(n), g(n), /(??), and k(n) for p=l and d € {1,2.3}. 25
3.3 Impulse responses h(n), g(n), /(n), and k(n) for p — 2 and c/ € {1,2,3}. 26
4.1 Impulse responses tx(n) for p € {0,1,2} 32
VII
LIST OF FIGURES
1.1 Discrete wavelet transform of two signals translated to each other. . . 2
1.2 Discrete wavelet transform and "algorithme a trous." 3
2.1 Spline functions p (x) and r) p (x) for p € {0,1,2,3,4} 10
2.2 Fourier transforms 3 p (ui) and rj p (u>) for p € {0,1,2,3,4} 13
2.3 Splines /5 p (x) and /3 p (w) for p € {0,1,2,3,4} 16
3.1 Filter bank implementation of a 1-D discrete dyadic wavelet transform. 21
3.2 Comparison of two discrete implementations using sinc(.r) as an input. 23
3.3 Wavelets r(.r) = d \j (x) for p € {0. 1,2} and d G {1.2.3} 27
4.1 Filter bank implementation of a 2-D discrete dyadic wavelet transform. 35
4.2 Filter bank implementation of a spline derivative-based transform. . . 48
6.1 Comparison of image fusion results 62
6.2 Comparison of image enhancement results 64
6.3 The magnitude of filter coefficients for G$(uJ x ,ujy) 65
vin
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
A MULTISCALE SPLINE DERIVATIVE-BASED TRANSFORM
FOR IMAGE FUSION AND ENHANCEMENT
By
Iztok Koren
December 1996
Chairman: Dr. Fred J. Taylor
Cochairman: Dr. Andrew F. Laine
Major Department: Electrical and Computer Engineering
Analyzing images along distinct scales is advantageous in many image pro-
cessing and computational vision tasks. Traditional wavelet transforms have become
a popular tool for multiscale image analysis and image compression suffer from arti-
facts, lack shift and rotation invariance, and may introduce aliasing and anisotropies.
We propose a highly redundant representation that eliminates these artifacts.
A new transform is constructed as an extension of the discrete dyadic wavelet
transform with a wavelet being the first derivative of a central B-spline. We begin
the construction in one dimension by providing an efficient initialization procedure
for the computation of the discrete dyadic wavelet transform and by extending the
transform to encompass higher order derivatives.
This modified one-dimensional discrete dyadic wavelet transform will serve
as a basis for deriving a steerable two-dimensional transform. Such a transform is
IX
advantageous in that it does not introduce any of the artifacts described above. How-
ever, exact reconstruction is problematic via implementation in the spatial domain.
We solve this problem by devising a spline-based approximation. The resulting algo-
rithm for computing a multiscale spline derivative-based transform is implemented
efficiently as a filter bank consisting of separable filters alone. For both finite and
infinite impulse response filters present in the filter bank, we specify implementa-
tion details of a realization that alleviates the boundary effects caused by the use of
circular convolution.
Finally, we use the derived transform for multiscale image fusion and enhance-
ment. We show empirically that the transform does not introduce artifacts commonly
reported for traditional wavelet-based implementations and provides more flexibility
for different fusion and enhancement strategies by enabling orientation analysis.
CHAPTER 1
INTRODUCTION
Analyzing images across multiple scales and resolution has become a power-
ful tool for solving compelling problems in computational vision, image processing,
and pattern recognition. Wavelet theory encompasses multiscale and multiresolution
representations, such as subband filtering [40], image pyramids [4], and scale space
filtering [48], into a unified mathematical framework. In the area of image processing,
there remain few research areas to which wavelet analysis has not been applied. For
example, problems in image compression, denoising, restoration, enhancement, reg-
istration, fusion, segmentation, and analysis, have all been approached with distinct
kinds of wavelet processing.
Though ubiquitous, wavelet analysis is not without problems of its own. Lack
of translation invariance, one of the major problems of the wavelet transform [10], is
in multiple dimensions accompanied with lack of rotation invariance.
1.1 Shortcomings of Traditional Methods of Wavelet Analysis
A wavelet transform in its most commonly used orthogonal or biorthogonal
forms is not translation and rotation-invariant. By translation-invariant transform,
we mean a transform that commutes with a translation operator. Since we will
deal primarily with discrete transforms in this work, we constrain the translation
parameter to integer multiples of a sampling period.
Lack of translation invariance of the discrete wavelet transform is illustrated
in Figure 1.1. Here, we can clearly see how a translation of the input signal by one
■
4 ill*
II
, !
i 4 in*
6
9 I r
(a)
(b)
Figure 1.1: (a) Original signal and (b) signal translated one sample to the left with its
discrete wavelet transform coefficients shown across dyadic scales 2 m , m G {1,2,3,4}.
sample results in a completely different set of transform coefficients (orthogonal
wavelet DAUB4 1 [10] was used in this experiment).
Noninvariance under translations of an orthogonal and biorthogonal wavelet
transform is due to lower sampling density at coarser scales. 2 A straightforward way
of dealing with this problem is to construct a redundant transform by using the same
sampling frequency for the input signal and all scales of the transform. A filter bank
implementation of such a transform, called "algorithme a trous" [15], is based upon
the fact that downsampling followed by filtering is equivalent to filtering with the
upsampled filter before the downsampling, as shown in Figure 1.2.
G(z) 2t
H(z) 2t
G(z) 2t
G(z)
G(z) 2t
H(z)
G(r)
H(z) 2t
(a)
G(z 4 )
H(r)
H(z')
(b)
Figure 1.2: Filter bank implementation for (a) a discrete wavelet transform and (b)
"algorithme a trous" decompositions for three levels of analysis.
Lack of rotation invariance is another shortcoming of traditional (i.e.. orthog-
onal and biorthogonal) wavelet techniques. In defining rotation invariance, we are a
bit less strict than with translation invariance. We do not require that the transform
commutes with a rotation operator here. Even in the case of a simple filtering, this
would limit us to circularly symmetric filters only. Our requirement for analysis is a
lr rhe number in DAUB4 refers to twice the order of the wavelet (i.e.. two in this case).
In practice, since analysis is performed over a finite range of scales, a discrete wavelet transform
is translation-invariant by translations determined by the coarsest scale (e.g., sixteen samples for
the analysis from Figure 1.1) [10].
transform that enables rotation-invariant processing. As an example of such a trans-
form, let us consider filtering with the first derivative of a two-dimensional Gaussian
probability density function in two directions, specifically, along x and along y-axis.
By linearly combining the results of filtering in these two directions, filtering with the
first derivative of a Gaussian in any direction can be computed. This fact was used
by Canny [6] for edge detection. A determined edge direction rotates as an input
image is rotated.
After choosing the fundamental properties of the transform, one must decide
upon the basis functions to be applied. For our studies, we selected basis functions
that well approximated derivatives of a Gaussian, because (1) the Gaussian proba-
bility density function is optimally concentrated in both time and frequency domain,
and thus suitable for time-frequency analysis, (2) higher order derivatives of a Gaus-
sian can be, similar to the first derivative, used for rotation-invariant processing [12],
and (3) the Gaussian function generates a causal (in a sense that a coarse scale de-
pends exclusively on the previous finer scale) scale space [2]. The last property makes
possible scale-space "tracking" of emergent features.
1.2 Thesis Overview
We wish to construct an invertible, translation and rotation-invariant two-
dimensional transform that decomposes an image using a set of basis functions which
are approximations to derivatives of the Gaussian probability density function.
To obtain a translation-invariant representation, we will use "algorithme a
trous," or more precisely, we will extend the discrete dyadic wavelet transform [28],
which uses "algorithme a trous." Derivatives of a Gaussian will be approximated by
wavelets which are derivatives of central B-spline functions. In Chapter 2, we review
the basics of B-spline processing necessary to derive transforms in later chapters.
Next, Chapter 3 examines a dyadic wavelet transform in one dimension. We
present extensions of the discrete dyadic wavelet transform to higher order derivatives
and to even order spline functions. To compute a discrete transform from a contin-
uous one, the discrete computation must be properly initialized. We devise a new
initialization procedure that is shown to be more accurate than the one suggested by
Mallat and Zhong [28]. After the derivation of the transform, we point out relevant
connections to scale-space filtering and reconstruction from edges.
In case of the first derivative of a Gaussian, a rot at ion- in variant transform
can be obtained by a tensor product extension of the one-dimensional transform
to two dimensions. For second order derivatives, a tensor product approach can
approximate Laplacian of a Gaussian, but cannot perform orientation analysis (due
to its isotropic nature). In Chapter 4, we first examine the above two cases, and
then develop a new multiscale spline derivative-based transform, which in addition
enables rotation-invariant processing for derivatives of orders higher than two.
Given that the developed transform is highly redundant, an efficient imple-
mentation is extremely advantageous. This issue is addressed in Chapter 5. When
filtering of a signal is performed by circular convolution, boundary effects may result.
A standard way of dealing with this problem is by mirror-extending the input signal
to the filter. We combine such an extension with symmetry/antisymmetry of the
filters, to achieve savings in both computation time and memory.
Finally, in Chapter 6, we present two sample applications: image fusion and
enhancement. We demonstrate that the new transform does not cause artifacts com-
monly associated with traditional wavelet techniques and is an attractive analytic
tool for designing novel fusion and enhancement strategies.
1.3 Notation
We use symbols TV, Z, and R for the sets of naturals, integers, and reals,
respectively.
L 2 {R) and L 2 (R ) denote the Hilbert spaces of measurable, square-integrable
functions f(x) and f(x,y), respectively.
The inner product of two functions f(x) £ L 2 (R) and g(x) £ L 2 (R) is given
by
(f(x),g(x)) = f°° f(x)g(x)dx.
J — oo
The norm of a function /(.r) £ L 2 (R) is defined as
11/11 = JJ~Jf(x)\*dx.
The convolution of functions /(x) £ L 2 {R) and </(x) £ L 2 (R) is computed as
/*</(-*•)= /" f(t)g(x-t)dt
J — 'OO
and the convolution of two functions f(x,y) £ L 2 {R 2 ) and g{x,y) £ L 2 (R 2 ) equals
f*g(x,y)= / f(tx,ty)g(x-t x ,y-t y )dt x dt y .
J—oo J—oo
The Fourier transform of a function /(.r) £ L 2 {R) is defined as
/(w)= / f(x)e-'»*dx,
J — oo
and the Fourier transform of a function f(x,y) £ L 2 {R 2 ) is equal to
/(u/„w y ) = r r f(x,y)e-^ + ^dxdy.
J —'OO J — oo
l 2 (Z) and 1 2 {Z 2 ) stand for the spaces of square-summable discrete signals f(n)
and f(n r ,riy). respectively.
The --transform of a discrete signal f(n) £ 1 2 {Z) is defined as
F(z)= £ f(n)z-\
n=—'Oo
The convolution of discrete signals f(n) £ P(Z) and g(n) £ / 2 (Z) is equal to
oo
f*g(n)= ^2 f( m )9(n-m),
m=— oc
and the convolution of discrete signals f(n x ,n y ) £ l 2 (Z 2 ) and g{n x ,n y ) £ 1 2 {Z 2 ) is
given by
oo oo
/ * g(n x ,n y )= Yl Yl f( m *> m y)g{ n r- m x , n y - m y ) .
m x =—oo m-u=—oo
The Fourier transform of a discrete signal f(n) £ l 2 (Z) is equal to the
transform evaluated on the unit circle
CO
F{u>)= £ /(»>-^ n ,
n= — oo
and the Fourier transform of a discrete signal f(n x ,n y ) £ l 2 (Z ) is defined as
FKw s )= £ £ /(n„n lr ) C -^- B - +< -» B »>.
For later use, we define the following functions:
1. the unit impulse function
r . . 1 for x =
I I) otherwise,
2. the unit step function
( 1 for x >
U{X) := for x < 0.
3. the rectangular function
4. the sine function
rect(x) := < n , ?
for x > i
• / v sin(7r.r)
smc(.r) := , and
1TX
5. the unit impulse sequence
£( \ J 1 ^° r ?? =
I otherwise,
where x £ jR and n £ Z .
CHAPTER 2
SIGNAL PROCESSING USING CENTRAL B-SPLINES
In this chapter we briefly review fundamentals of spline processing needed for
derivations in Chapters 3 and 4.
2.1 Central B-Splines: Definition and Properties
Given real numbers — oo < x < x\ < ,r 2 < . . . < x m < x m+ i < oo, a function
on the interval [.ro,.r m +i] is called a spline function of order p with the knot (i.e.,
grid point) sequence Xi, x-2. ... .r m . if it is (1) a polynomial of degree p or less in each
interval [a?,-, x i+ i]. i = 0, 1, . . . m, and (2) continuous in its derivatives up to the order
p-1 on the interval [.r ,.r m +i] (i.e., C p ~ 1 [x ,Xm+i})-
Here, we will concentrate primarily on basis splines (B-splines), or more pre-
cisely, central B-splines having knots at i £ Z for p odd and at i + | for p even [34].
Central B-splines of order p (with p+1 knots) are defined as
Figure 2.1 shows 3 p (x) and their Fourier transforms (5 p {u)) for p £ {0, 1, 2, 3, 4}.
A family of functions {/3 p (x — ???)} % forms a basis of S p , a space of order
p spline functions with knots at i for p odd and at / + | for p even (?' G Z) [34. 35].
Except for p = 0, the basis functions {/? p (x — m)} are not orthogonal.
Let us list some properties of functions j3 p [x) [3, 34]:
1. fip{x) are nonnegative functions with a support of length p+1,
p-\-l times
/ * V
2. /# p (.r) = fio * /?o * • • • * A)(#) 5 where "*" denotes the convolution operator, or,
equivalently, in the Fourier domain:
3. P P {x) = I ((*? + x) & P - X {x + \) + (*!-*) ^(x - §)),
*-¥ = M j+ !)-M*-!)-
Another interesting property of B-splines is the fact that they converge to a
Gaussian as their order tends to infinity. Unser et al. [44] derived the Gaussian
approximation
i -4
where a v — J 2 ^-- Ratio n7 ,n , for example, is already well below 0.2%.
v i^ lift* 1 1
Denoting by S\ a spline function space spanned by {fi p (x — m)} for p odd
and by {/3 p (x — \ — rn)} for p even (subscript in Sj refers to the fact that spline
functions have knots at integers), spline spaces form a nested sequence . . . C St, C
. . . C S j C Sj C S^-i C . . . C S^-i C .... By orthogonalizing this basis functions
a multiresolution analysis of L 2 (R), from which the Battle- Lemarie wavelet bases
stem, can be built [10]. Here, we will not pursue constructions of orthogonal, semi-
orthogonal, or biorthogonal spline-wavelets any further — reader looking for a detailed
treatment of this subject may find a good starting point in books [8, 9].
2.2 B-Spline Signal Interpolations
Unser et al. developed a fast digital filtering scheme for B-spline signal pro-
cessing [43]. They defined a discrete B-spline of order p and expansion factor (spacing
between knots) m as
b p ,m(n) := ,8 P l—j , n,m € Z.
10
p-n
-3 -2 -1
I a |
(b)
Figure 2.1: Spline functions (a) j3 p (x) and (b) t} p (x) for p <E {0, 1,2,3,4}.
II
Henceforth, if the distance between knots equals one, we will write b p (n) instead of
Interpolation of a discrete signal s(n) G l 2 (Z) by s p (x) G S p using central
B-splines
s(n) = s p {x)
J2 c(i)/3 p (x-i)
5
x=n
(2.1)
can now be written as a convolution sum
s(n) = c*b p (n). (2.2)
If s(n) are samples of a function s(x) bandlimited to [— k. k] (i.e., the support
of its Fourier transform S(uj) is in [— 7r,7r]), it can be shown that s p {x) — ► s(x) as
p — > oo [1, 36].
In [13]. they refer to a linear operator by which B-spline coefficients c(n) can
be obtained from samples s(n) as a "direct B-spline transform." Equation (2.2)
therefore represents "indirect B-spline transform'' of a sequence {c(n)}.
After taking the --transform of (2.2), the direct B-spline filters are found to be
^p 1 ! 77 ) = Z~ 1 {[Bp( z )]~ 1 }- Since B-spline functions j3 p (x) are compactly supported,
indirect B-spline filters b p (n) are finite impulse response (FIR) filters, while direct
B-spline filters ^(n) are infinite impulse response (IIR) filters. Aldroubi et al. [1]
showed that IIR filters b~ l {n) are stable (i.e., the region of convergence of B~ l {z)
includes the unit circle [30] ) for any order p. Note that both indirect and direct B-
spline filters are symmetric, which follows from the fact that central B-splines fi p {x)
are symmetric.
Table 2.1 shows the --transforms of direct B-spline filters for the first ten
orders. We postpone the discussion on implementation details of B-spline filters until
Chapter 5.
L2
Table 2.1: Transfer functions of direct B-spline filters for orders from to 9.
vo
8
6
z+4+z -1
384
~? +76J+230+76J- 1 +z~ 2
120
j 2 +26-+66+26;- 1 +.— 2
46080
'4-722; 2 + 10543j4-23548+10.543j- 1 +T22j- 2 +j- 3
5040
j 3 + 120; 2 + 1191 i+2416+l 191 z~* +120z~ 2 +z- 3
10321920
: 4 +6552j 3 +331612- 2 +24852883+4675014+2485288;- 1 +331612z- 2 +6552^- 3 + ;- 4
3(12880
: 4 +502- 3 + 14608s 2 +88234z+156190+88234z- 1 +14608;- 2 +502:- 3 + -- 4
Instead of using B-spline interpolation as given by (2.1) it is sometimes con-
venient to express the interpolating function s p (x) in terms of discrete samples s(n)
s p( x ) = Y, s ( l )Vpi x -i),
(2.3)
where r/ p (.r ) is the cardinal spline of order p. In the frequency domain, cardinal splines
converge to an ideal lowpass filter with cutoff frequency tv (i.e., r] p (x) — ► sinc(.r)) as p
tends to infinity [1. 45], which establishes the asymptotic equivalence with Shannon's
sampling theorem [37].
Using (2.1) and (2.2) with (2.3) cardinal splines can be related to B-splines:
V P (x)= E b;\i)P P (x-i). (2.4)
i=—oo
Cardinal splines T) p (x) and rj p {w) for p G {0, 1,2,3,4} are shown in Figure 2.2.
13
-2s 2s
-2» 2s
-8s -6s ^ -2s 2s 4* 6a 8n
-6s -4s -2s 2s 4™ 6* 8s
(a)
-6s -4* -21 2s 4: 6x 8*
-8s -6s -4s -2 7 2* 4s 6s 6*
-8s -6s -4* -2* 2s 4s 5j
-6s -6s -4s -2s 2s 4s 6s 8s
-8s -6s -4s -2s 2. 4s 6s 8s
(b)
Figure 2.2: Fourier transforms of spline functions (a) 3 p {lj) and (b) rjJu) for p £
{0,1,2,3,4}.
II
2.3 B-Spline Signal Approximations
Central B-splines are also simple to use when the goal is function approxima-
tion. Least-squares B-spline approximation of s(x) £ L 2 {R) is achieved by computing
the orthogonal projection of this function onto S p . We have
~s(x) = s p (x) = f] d(i)/3 p (x-i) (2.5)
i=—oo
with
d(i) = (s(x)J p {x-i)), (2.6)
where
°°
&(*) = £ a&iCO&C* -
!= — OO
O
is the dual spline of order p [45]. Spline functions /3 p (x) and /? p (:r) satisfy the biorthog-
onality condition
o
((3 p (x - m), i3 p (x - n)) = 6(m - n), m, n e Z.
(Note that, since both /3 p (x) and ft p (x) form a basis of 5 P , they can be interchanged
in (2.5) and (2.6).)
O o
Figure 2.3 shows functions f3 p (x) and their Fourier transforms 3 p (uj) for p 6
{0,1,2,3,4}.
An interesting alternative to the minimum Z 2 -norm (i.e., least-squares) ap-
proximation of a signal is obtained by computing oblique instead of orthogonal pro-
jection of the signal onto the spline function space. Unser and Aldroubi proposed an
independent specification of the sampling and approximation spaces [42]: a linear op-
erator maps coefficients of the input signal expansion over sampling space basis into
the coefficients of the approximation space basis expansion from which the signal's
projection onto the approximation space is recovered. Constraining the entire system
to be linear, shift-invariant for integer translations, and consistent (i.e., the system
acts as an identity operator for functions that belong to the approximation space),
l.-J
the obtained solution for the signal approximation is the projection of a signal onto
the approximation space perpendicular to the sampling space [42]. Analogously to
(2.5) and (2.6) this projection can be expressed as
sr(x)= f] a{i)0 r (x-i) (2.7)
i=— oo
with
a(i) = q£ *(s(x),p 3 (x-i)),
where q^} is the convolution inverse of the cross-correlation sequence
?«(») = {P$(x - i),Pr(x)) = b s+r+1 {i).
When the sampling space S s and the approximation space S r are identical
(i.e., 5 = r), an orthogonal projection given by (2.5) and (2.6) results. (Note that
the described oblique projection is not restricted to spline function spaces — the only
requirement is that both sampling space basis and approximation space basis are
Riesz bases of the corresponding function spaces [42].)
Signal approximation (2.7) is particularly attractive in situations where the
sampling space is given a priori (e.g., by the impulse response of the acquisition device
[42]) or when such an approximation is close to the optimal least-squares solution but
simpler to implement than the orthogonal projection (e.g., [47]).
I(i
-15 -to -5
-15 -tO -5
S 10 IS
5 10
5 tO 15
10 15
0-1
s
■
-8* -6« -4* -2i ir
P-2
. _J
,
-At 2* 4» &> Si
-<*I -2* 2i 4i 6i b
5 10 15
(a)
o
Figure 2.3: Spline functions (a) /3 p (x) and (b) their respective Fourier transforms
4(«)forp€ {0,1,2,3,4}.
CHAPTER 3
TRANSFORM IN ONE DIMENSION
In this chapter, we will augment the one-dimensional discrete dyadic wavelet
transform [28] to obtain a strong foundation for derivations of two-dimensional trans-
forms in Chapter 4.
3.1 1-D Discrete Dyadic Wavelet Transform Revisited
A discrete wavelet transform is obtained from a continuous representation
by discretizing dilation and translation parameters such that the resulting set of
wavelets constitutes a frame. The dilation parameter is typically discretized by an
exponential sampling with a fixed dilation step and the translation parameter by
integer multiples of a fixed step [10]. Unfortunately, the resulting transform is variant
under translations, a property which makes it less attractive for image analysis (cf.
Section 1.1).
As we have already mentioned in Chapter 1, sampling the translation param-
eter with the same sampling period as the input function to the transform results in
a translation-invariant, but redundant representation. The dyadic wavelet transform
proposed by Mallat and Zhong [28] is one such representation. Let us begin with a
brief review of properties of the dyadic wavelet transform as described in [28], but
included here for completeness.
The dyadic wavelet transform of a function s(x) € L 2 (R) is defined as a
sequence of functions
nM^)} rae z- (3-1)
IT
IN
where W m s(x) = s * i^ m (x), and m (.T) = 2 -m ^(2 -m x) is a wavelet tp(x) expanded
by a dilation parameter (or scale) 2 m . Note the use of convolution instead of an inner
product.
To ensure coverage of the frequency axis the requirement on the Fourier trans-
form of V'm(-c) is the existence of A\ > and B\ < oo such that
oo
A, < £ |0(2 m cu)| 2 < B x
m=— oo
is satisfied almost everywhere. The constraint on the Fourier transform of the
(nonunique) reconstructing function \(x) is
m=— oo
A function s(x) can then be completely reconstructed from its dyadic wavelet trans-
form using the identity
oo
S ( X ) = Yl W ^ S * Xm(x),
m=— oo
where \ m {x) = 2- m \(2- m .r).
In numerical applications, processing is performed on discrete rather than
continuous functions. When the function to be transformed is in the discrete form,
the scale 2 m can no longer vary over all m £ Z. Finite sampling rate prohibits the
scale from being arbitrarily small, while computational resources restrict the use of
an arbitrarily large scale. Let the finest scale be normalized to 1 and the coarsest
scale set to 2 .
The smoothing of a function s(x) G L 2 {R) is defined as
S m s(x) = s*<f> m (x),
where <f> m (x) = 2~ m (p(2~ m x) with m G Z, and <j>(x) is a smoothing function (i.e., its
integral is equal to 1 and <f>(x) -+ as Ixl — > oo).
1!)
In [28], a real smoothing function (f>(x) was selected, whose Fourier transform
satisfied
|^)| 2 = f)^(2 m w)x(2 m o;). (3.2)
m=l
In addition, it was shown that any discrete function of finite energy s(n) £ 1 2 {Z)
can be written as the uniform sampling of some function smoothed at scale 1, i.e.,
s(n) = Sof(n), where f(x) G L 2 (R) is not unique. Thus, the discrete dyadic wavelet
transform of s(n) for any coarse scale 2 M was defined as a sequence of discrete func-
tions
{S M f(n + a), {W m f(n + s)} m£[hM] } neZ ,
where s is a ir{x) dependent sampling shift.
The above initialization s(n) - S f(n) is rather standard in the discrete
wavelet transform computation [10]. although it yields correct results (i.e., the dis-
crete wavelet transform is equal to the samples of its continuous counterpart) only
when s(n) = Sqs(ti). Here, we will concentrate on wavelets which are derivatives of
spline functions and this will lead us to a simple initialization procedure [46] that
alleviates the above problem.
For a certain choice of wavelets the discrete dyadic wavelet transform can
be implemented within a fast hierarchical digital filtering scheme. Next, we shall
summarize the relations between filters, wavelets, and smoothing functions.
First, let us introduce a real smoothing function <p(x) such that (3.2) can be
rewritten as 1
oo
#(w)0M= £^(2"V) x (2 m a;), (3.3)
m=0
and let us set 6{x) = /3 p (x) (i.e., we restrict ourselves to wavelets which are spline
functions).
^ote that the sum index determines the range of scales of the discrete transform: using (3.2)
we have ip(2u>) and \(2u) at the finest scale of the transform, while for (3.3) we get ip(u) and \(u).
20
Computing (3.3) for the finest two scales shows that
0(w) x(u) = A(«)^(«) - &(2w) 0(2w). (3.4)
/5 p (2w) can be related to /i p (u;) by expressing /i p (2u>) as (cf. Proposition 1 of [46])
1 / sin(^) \ / sin UJ
p (2w)
2P+1 \™(l).
and using ELo^'"^ = !i ^lf 1 e4 ^ ):
4,(2w)=(oob(|)) P+1 A(«). (3.5)
(Note that a relation similar to (3.5) can be derived for integer scales provided that
the dilation parameter and the order p are not both even [46].)
Let F(m) be a digital filter frequency response and let F s (uj) = e JUJS F(io).
If we choose
0(w) = G_ # (w)4p(«), (3.6)
tp(2u)=L,(u)<p(u) t (3.7)
x(u) = K.(u)<p(w), (3.8)
and
H(u:) = e^(cos(f)Y'\ (3 " 9)
where s G {0. |} is a filter dependent sampling shift needed for g(n). /(/?). k(n), and
h(n) to be FIR filters, and insert Equations (3.5)-(3.9) into (3.4), we observe the
relation between frequency responses of the filters
G(u)K(u) + H(u)L(u) = 1. (3.10)
Similar to orthogonal and biorthogonal discrete wavelet transforms, the dis-
crete dyadic wavelet transform can be implemented within a hierarchical filtering
21
scheme. To derive such a digital filtering scheme, let us assume that s(iv) from (3.1 )
is bandlimited to [— 7r,7r]. Using Shannon's sampling theorem [37] and (3.6) in the
definition of the dyadic wavelet transform (3.1) with m = shows
II u
/oo °° °°
J2 s ( i )sinc( t - i ) ]T g_ s ( m ) 3 P ( x - t - m ) dt .
-oo • _, _,_
By making use of the fact that the cardinal spline functions tend to the sine
function as their order r approaches infinity, and employing (2.4) we can write
Wo»H * S{u>)B?{u)P r (u>)$ p {a,)G-.{u,),
or, by using (3.5) and (3.9),
m — 1
^{VU m .(.r)| x= J~.SV)^- 1 ^)i? P+r+ i^)G- 5 (2 m a;)n^- s (2 ! -'). (3.11)
;=o
Equation (3.11) entirely specifies the discrete dyadic wavelet transform de-
composition, while the reconstruction follows from (3.4)-(3.9). Three levels of a filter
bank implementation are shown in Figure 3.1. (Note that the initialization is the
same as the one proposed in [46] and that except for the prefiltering and postfilter-
ing, this scheme is implementing "algorithme a trous.") Noninteger shifts at scale 1
for filters with s = | are rounded to the nearest integer.
— ^o» BpWco)
G>)
- H«o)
K,(a,)
Q(2<0)
H(2<o)
G.(4<o)
HJ4<:,)
^(4(0)
M4co)
1^(20.)
L s (2<o)
M»)
^,«o); 8,(0.)
Figure 3.1: Filter bank implementation of a one-dimensional discrete dyadic wavelet
transform decomposition (left) and reconstruction (right) for three levels of analysis.
We will now perform a simple experiment which will illustrate the difference
between the implementation of the discrete dyadic wavelet transform as originally
22
proposed in [28] (i.e., without prefiltering and postfiltering) and the one from Figure
3.1.
Let s(x) = sinc(.r), p — 2, and g(n) = 28 (n + 1) — 28(n) (this particular choice
for p and g(n) results in the same wavelet as was used by Mallat and collaborators
[27, 28]). The dyadic wavelet transform of s(x) at a scale 2 m (3.1) in the frequency
domain is then
C*(w) = G. s {2 m u) fa(2 m u) rect (^) . (3.12)
The Fourier transform of the discrete dyadic wavelet transform of s(n) = 8(n) at a
scale 2 m using spline based initialization follows from (3.11)
m-l
F{W m s(n)} = B;\«:)B r+3 ^)G- s (2 m u,<) Y[H_ S (2'^). (3.13)
i=0
while the one using the algorithm from [28] equals
m — 1
F{W m s{n)} = G-,(2 m u>) J] H..(Tu). (3.14)
i=0
In Figure 3.2 a comparison of magnitudes of (3.13) and (3.14) versus (3.12) is
shown: in Figure 3.2(a) magnitudes of (3.12) (solid) and (3.14) (dashed) are plotted
for m € {0. 1,2.3}, while the dashed curves in 3.2(b) represent magnitudes of (3.13)
with r = 5.
By choosing the appropriate order r, (3.13) can approximate (3.12) in the
interval [ — 7r, tt] arbitrarily good, while originally proposed (3.14) has troubles at finer
scales. Mallat and Zhong [28] noticed that there was a problem with their discrete
transform computation, and introduced a set of constants associated with the discrete
transform coefficients at dyadic scales. They chose the values of constants such that
the transform coefficient modulus maxima remained constant over all dyadic scales for
a step edge input signal. In relation to Figure 3.2(a) this is equivalent to multiplying
T{W m s(n)} by a distinct constant for each m. Clearly, this can improve over the
situation depicted by Figure 3.2(a), but can still not compete with the described
spline based initialization.
23
(a)
(b)
Figure 3.2: (a) Fourier transform magnitudes of the dyadic wavelet transform of
s(x) = sinc(x) (solid) and the originally proposed discrete dyadic wavelet transform
[28] of s(n) = S(n) (dashed), (b) Fourier transform magnitudes of the dyadic wavelet
transform of s(x) (solid) and the discrete dyadic wavelet transform using quintic
splines for interpolation of s(n) (dashed).
24
Next, we will choose filters in the filter bank implementation of the discrete
dyadic wavelet transform. As already mentioned, we are interested in wavelets which
are derivatives of spline functions. G(u>) in (3.6) is therefore the Fourier transform of
the difference operator centered around —3 (cf. Property 4 in Section 2):
G(w) = e*" 2; sin -J . (3.15)
dmod2
where d is the order of the derivative and the sampling shift for this filter is „
Since H(uj) was already given by (3.9), the remaining two filters to be deter-
mined are L(u>) and K(ui). Both of them are (as is true for <p(x) and \{x)) nonunique.
If we choose L(u>) such that we can express K[u) in terms of a finite geometric
series having the smallest number of elements for an arbitrary p, we get
L( U ) = e-^ y. (-i)- +l [ l ^ j \ ^ cos y j (3 16)
and
i d+1 <
/ //.A \ dmod'2 / P / /, ,\ \ 2ra\ L 2
' ,-Jujs ■ I w *
where [x\ denotes the largest integer smaller than x, the sampling shift for L{u)) is
the same as the one for H{u) (i.e., a= ^)^ od2 ) 1 and the sampling shift for K"(w) is
the same as the one for G(u)).
Note that Equations (3.9) and (3.15)-(3.17) work fine from the mathematical
point of view, but in practice the reconstruction may become cumbersome when both
p and d are large (the lengths of impulse responses ft(n), g(n), l(n), and k(n) are p+2,
d+1, (p+l)(d-(d+l)mod2)+l, and pd+(p+l)(dmod2)+l, respectively, while for the
frequency responses of the decomposition filters we observe that linip^oo \H- S {uj)\ =
8 u (u + 2nir) and Yim d ->^(2j)- d |GL,(w)| = S u (u + {2n+l)ir) with n € Z).
It is also worth noting that if(w) is a lowpass filter when p is even (i.e., the
reconstruction function \{x) is a wavelet only for p odd).
25
Tables 3.1, 3.2, and 3.3 list impulse responses of the four filters for p 6 {0, 1, 2}
and d 6 {1, 2, 3}, while Figure 3.3 shows wavelets 0(x) = ' j*d f° r the same values
of p and d. Wavelets from this family have a support of length d+p+1, regularity
order p, and are either symmetric or antisymmetric.
Table 3.1: Impulse responses h(n), g(n), l(n), and k(n
n
h(n)
d= 1
d = 2
d = 3
g(n)
l(n)
k(n)
9(n)
l(n)
k(n)
g(»)
l(n)
k(n)
_2
i
-1
0.5
1
1
-3
-0.125
0.5
-1
0.5
-0.25
-2
0.5
-0.25
3
0.625
0.0625
1
0.5
0.25
1
0.5
-1
0.625
-0.0625
2
-0.125
for p = and d € {1,2, 3}.
Table 3.2: Impulse responses h(n), g{n), /(??), and k(n) for p=l and d £ {1,2,3}.
n
h(n)
l(n)
d = 1
d = 2
i/(»)
k(n)
9(n)
k(n)
-1
0.25
1
-0.0625
1
-0.0625
0.5
-1
-0.3125
-2
-0.375
1
0.25
0.3125
1
-0.0625
2
0.0625
n
J = 3
fe(n)
g(n)
l(n)
fe(n)
-3
-0.015625
-2
1
-0.09375
0.00390625
-1
0.25
-3
0.265625
0.04296875
0.5
3
0.6875
0.1015625
1
0.25
-1
0.265625
-0.1015625
2
-0.09375
-0.04296875
3
-0.015625
-0.00390625
As already discussed, wavelets with p = 2 and d=\ from a family of wavelets
with p even and d=\ were used in [27, 28], whereas filters with p= 1 and d = 2 from a
family of filters with p odd and d = 2 were employed by Laine and collaborators [1 1 , 22,
23]. Here described transform puts no restrictions on the choice of p or d whatsoever,
and uses a better initialization procedure than the one originally proposed in [28].
2(i
Table 3.3: Impulse responses h(n), g(n), /(»), and k(n) for p — 2 and d (E {1,2,3}
7?
fc(n)
d =
1
d =
2
9(n)
l(n)
k(n)
9(n)
'(»)
fc(n)
-2
0.125
-0.015625
-0.015625
-1
0.375
1
0.125
-0.109375
1
0.125
-0.125
0.375
-1
0.375
-0.34375
-2
0.375
-0.46875
1
0.125
0.375
0.34375
1
0.375
-0.125
2
0.125
0.109375
0.125
-0.015625
3
0.015625
n
d = 3
&(n)
g{n)
l(n)
fc(n)
-4
-0.001953125
0.000244140625
-3
-0.017578125
0.003662109375
-2
0.125
1
-0.0703125
0.0263671875
-1
0.375
-3
0.085937
0.0908203125
0.375
3
0.50390625
0.13037109375
1
0.125
-1
0.50390625
-0.13037109375
2
0.0859375
-0.0908203125
3
-0.0703125
-0.0263671875
4
-0.017578125
-0.003662109375
5
-0.001953125
-0.000244140625
27
-0.2
-0 6
---71-
— ,
-
■
]/ ~~\
\\ |
'
-
/
1 '
.
-
<'
\\
b
-3
-2
-1
2
3
(a)
1
-
--jr-
,__ .
-
05
A
\\
\l
V
J
\\
i
V
0.5
-
-1
V 7
\\ //
-
1 5
-
-2
-
Figure 3.3: (a) Wavelets (a) 0(x) = ^ M , (b) 0(ar) = ^ff^ , and (c) 0(*)
d3/? g; 3 3(j) for p = (dashdotted), p=l (solid), and p = 2 (dashed).
29
the zero-crossings of W m s(x) for d = 2, a scale-space image of the representation sim-
ilar to [48] can be obtained. The only differences are that curves in the scale-space
plane are computed for dyadic scales only and that i) m is a close approximation of a
Gaussian instead of a true Gaussian.
CHAPTER 4
TRANSFORMS IN TWO DIMENSIONS
In this chapter, we will extend the one-dimensional transform from Chapter
3 to two dimensions 1 and derive several two-dimensional transforms with translation
and at least approximately rotation-invariant decomposition.
4.1 2-D Discrete Dyadic Wavelet Transform Revisited
The dyadic wavelet transform of a function s(x,y) G L 2 (R 2 ) is defined as a
set of functions [28]
{Wls(x,y),Wls(x,y)} m€Z , (4.1)
where W^s(x,y) = s * ^ m (x,y) for i = {1,2} and ^ m (x,y) = 2- 2m t l {2- m x,2- m y)
are wavelets ^'(x,y) expanded by a dilation parameter 2 m .
To ensure coverage of the frequency space there must exist an A 2 > and
B 2 < oc such that
A 2 < £) j:iV'*'(2 m ^,2 m u; v )| 2 < B 2
m = — oo i=\
is satisfied almost everywhere. If (nonunique) functions \ l {x.y). \ 2 (x,y) are chosen
such that their Fourier transforms satisfy
£ ^^(2^,2^)^(2^,2-^) = 1.
l = — oo z'=l
the function s(x,y) may be reconstructed from its dyadic wavelet transform by
4^y)= f; s:wi«*xi.(*,y).
7n=— oo i=l
where \j n (.r,y) = 2~ 2m x i (2~ m x,2- m y).
1 For extensions to higher dimensions, please refer to [18, 19].
30
31
However, when processing discrete functions the scale 2 m may no longer vary
over all m € Z . Let the finest scale be normalized to 1 and the coarsest scale set to
be 2 M . Let us, similar to [28], introduce a real smoothing function (f>(x,y) such that
its Fourier transform satisfies
|#c*,<*)| a = £ E^(2 m ^,,2"H) \'C2 m ^2 m ^ y ). (4.2)
m =0 /=1
Here, as in one dimension, a finite energy discrete function (s(n x ,n y ) £ l 2 (Z )) can
be written as the uniform sampling of some function smoothed at scale 1: s(n x ,n y )
= S f(n x ,n y ), where f(x,y) G L 2 {R 2 ) is not unique, and S m f(x,y) = /* <f> m (x,y).
This led Mallat and Zhong [28] to a two-dimensional analog of the one-dimensional
definition of the discrete dyadic wavelet transform: 2
{ S M -if{n x +s, n y +s),{W^f(n x +s,n y +s),W^f(n x +s,ny+s)} me [o t M-i]} {nxny)&Z 2-
We will use, as in Section 3.1, a spline-based initialization procedure.
To implement a multidimensional discrete dyadic wavelet transform within a
fast hierarchical digital filtering scheme, the wavelets were chosen to be separable
products of one-dimensional functions:
1 (*,y) = ^(ar)^(y), (4.3)
4>\x,y) = V(y) </>(x), (4.4)
where d(.r) and 0(.r) were chosen as described in Section 3.1 (i.e.. d(x) = /3 p (x) and
^>(u>) specified by (3.6)).
From (4.3), (4.4), and (3.6), we may write
frfatUy) = G- s {u! r )l3 p {u T )l] p {u;y), (4.5).
j> 2 (u x ,uj y ) = G- S {u y )f3 p {u x )f3 p {uj y ), (4.6)
2 As in Section 3.1, we put the finest scale of the transform at »7) = 0.
32
where G(w) is given by (3.15) for d G {1,2}. Choosing
A A
X (u x ,u> y ) = K a (u x )Ti(wy)P P (u) a .)P P (u) y ),
X 2 (U X ,U}y) = K a (LJy)Tl(U X )P p (u) X )0 p (u}y),
(4-7)
(4.8)
where K(u>) and Ti(uo) are digital filter frequency responses, we may compute (4.2)
for the finest two scales by
X> f (2-v,2^) x\2u t ,2u y ) = \Uu x ,ui v )\ 2 - |#(2w„2w,)| 3 . (4.9)
, = i
Inserting the terms defined by (4.5), (4.6), (4.7), (4.8), (3.5), and (3.9) with <j>{ut x ,w v ) =
^ p {u} x )^ p (ijJy) into (4.9) results in
K{uj x )G{ui x )T x {u v ) + KfaWMTrfe) + \H(u x )\ 2 \H(u y )\ 2 = 1. (4.10)
Equation (4.10) represents a relation between the frequency responses of the digital
filters used to implement a multidimensional discrete dyadic wavelet transform and
is a multidimensional analog to (3.10).
Solving (4.10) for 2i(u;) by substituting K(uj)G(u)) from (3.10) yields the
closed formula [28]
T^)= [ -(l + \H(^\ 2 ) (4.11)
In Table 4.1 we provide the filter coefficients for Ti(uj) from (4.11) for p £ {0, 1,2}.
All other filters from (4.10) were already specified in Section 3.1.
Table 4.1
: Im
Dulse responses t
i(n) for p e
n
p=0
p=l
p=2
-3
0.0078125
-2
0.03125
0.046875
-1
0.125
0.125
0.1171875
0.75
0.6875
0.65625
1
0.125
0.125
0.1171875
2
0.03125
0.046875
•J
0.0078125
33
As in the one-dimensional case, a two-dimensional discrete dyadic wavelet
transform can be implemented as a fast hierarchical filtering scheme. To derive
such an implementation, we, similar to the one-dimensional case from Section 3.1,
use the definition of the two-dimensional dyadic wavelet transform (4.1) and require
s(uj x ,u} y ) = for \uj x \ > 7r or \u> y \ > k. Using Shannon's sampling theorem in two
dimensions [16], (4.5), (4.6), and m = 0, we get
oo oo
/oo /-co -Q -"-'
J2 Yl s(i x ,iy)smc(t x -i x )smc(t y -i y y
-CO J — CO ■
l x =-rx., l y = -,yz>
X
• Y, 9-s{m)(3 p (x -t x - m)f3 p {y - t y ) dt x dt y
n
and
/CO /-oo ' — ' -'-'
Y. H s{i x ,i y )smc(t x -i x )smc(ty
■° oJ -° i x =~ooi v =- o
oo
•Pp(x - t x ) J2 9-s(m)P P (y -t y - m) dt x dt y .
m — — co
oo oo
111 = — CO
As in one dimension, we make use of the fact that the cardinal spline functions
converge to the sine function as their order r tends to infinity, and write
W^s(uj x ,u y ) ~ S(u x1 u y )B~ 1 (u} x )B~ 1 (u}y)^ p+r+1 {u} x )^ p+T+1 (u; y )G- a (ut x )
an(
W$s(u} x ,u y ) ~ S(u t ,u y )B r 1 (u x )B r 1 (u}y)$ p+r+ i(u x )^ p+r+1 (u} y )G- s (u; y ),
or for m £ [0, M) and discrete signal processing
F{W l m s(x,y) } ~ S(u x ,ujy)B;\u; x )B^(u;y)B p+r+1 (u x y
x=n x ,y=n v
m—1
■B p+r+1 (u y )G^(2 m uj x ) J] H_ s (?u: x )H_ s (?uJ y ). (4.12)
t'=0
and
F{W 2 m s{x,y) _ _ }~S(u x ,u y )B;\u x )B; 1 (u y )B p+r+1 (u, x y
771 — 1
■B p+r+1 (u y )G. s (2 m i, y )l[H^C2 i ^.)H_ s C2^ y ). (4.13)
i=0
Equations (4.12) and (4.13) describe the decomposition part of the filter bank
implementation of a two-dimensional discrete dyadic wavelet transform. The recon-
struction part can be obtained from (4.7)-(4.9) with (3.5) and (3.9). The entire filter
bank implementation of the transform is shown in Figure 4.1. Except for the pre-
filtering and postfiltering, we readily recognize the implementation proposed in [28].
Using the fact that a wavelet V'( x ) ^ s equal to a first (d= 1) or a second {d = 2)
derivative of a spline function /? p+ d(x), (4.3) and (4.4) may be rewritten as
V[x,y) = — -Q^i — , i,d€{l,2],
where
i) i (x.y) = (3 p+d (x)My)
and
d 2 (x,y) = P P (x)p p+d (y).
Let us denote W m s(x,y) = (W^s(x, y), W^s(x. y)). V = (£, £), A = V 2 =
(|^ + Jt-), and assume that d'(x,y, ) can be approximated by O(.v.y) for both / £
{1.2}.
For c/ = 1 it then follows that
W m s(x,y) = 2 m V{s*# m ){x,y). (4.14)
Thus for c/ = 2 we can write
2
'£Wi nS (x,y) = 2 2m A(s*# m )(x,y). (4.15)
1=1
With ${x, y) being a Gaussian, finding local extrema of (4.14) in the direction
of gradient V corresponds to the filtering stage of a Canny edge detector [6], and
finding zero-crossings of (4.15) corresponds to the filtering carried out with a Marr-
Hildreth edge detector (Laplacian of Gaussian) [29]. (Note that both edge detectors
Bftojlgffly) B p+r+1 (m x ) B^ !(cp y )
G-s(co x )
G s (co v )
H s (co x ) H^)
(a)
G,(2co x )
G s (2 My )
35
H s (2co x )H s (2o) y )
K s (co x ) T^coJ
1,(0),) K s (o) )
Ks(2co x ) 1,(2(0,)
T 1 (2co x )K s (2 My )
'H
^(2^)^(2^)
H -sK) H -*(co y )
(b)
B p 1 +r+l (o) x )B p 1 +r fl (cQ y ) -- B r K) B r (co y )
Figure 4.1: Filter bank implementation of a two-dimensional discrete dyadic wavelet
transform (a) decomposition and (b) reconstruction for two levels of analysis. H* s (u)
denotes the complex conjugate of H- s (uj).
36
involve postprocessing). Edge detection based on finding local extrema of W m s(x.y)
or zero-crossings of £2, =1 W^ n s(x,y) is therefore an approximation to the Canny or
the Marr-Hildreth edge detector over a range of dyadic scales. The differences stem
from the fact that r)(x.y) is neither a Gaussian nor is i?*(x,y) equal to i)(x,y).
4.2 Steerable Functions
When extending the transform from Chapter 3 to two dimensions, one of the
first questions that conic to mind is how to deal with the fact that derivatives can
be defined in any direction of the plane. 3 In case of a first derivative of a Gaussian,
one can simply compute first derivatives of a Gaussian in x and y directions and
then determine the derivative in any direction from these two directional derivatives
[6]. For higher order derivatives of a Gaussian, Freeman and Adelson [12] showed
that order+1 directional operators are needed for synthesizing the operator at any
orientation. In fact, functions with the property of expressing their arbitrary rotations
as linear combinations of some functions are not limited to derivatives of a Gaussian.
Below, we briefly restate some of the results from [12].
A two-dimensional function is called "steerable" if its rotations generate a
finite dimensional space. Rotations of a steerable function f(i\0) can therefore be
expressed as
/M-0 o ) = £ci(0o)/,(r,0), (4.16)
where 6 denotes an arbitrary angle, {cj(0 o )} stands for a set of interpolating func-
tions, {/,(r, 0)} is a set of basis functions, and r = y/x 2 + y 2 and 6 = arg(.r,y) are
polar radius and angle, respectively.
If /(r, 6) represents a filter kernel, the result of filtering with a rotated fil-
ter f{r\6 - 6 ) can be computed simply by {c;(0 o )} weighted linear combination of
outputs from basis filters {f t (r,8)}. Only the outputs from basis filters need to be
3 Second derivatives of central B-splines can be used, as we saw in Section 4.1, to approximate
Laplacian of a Gaussian for approximately rotation-invariant, although not directional, processing.
M
precomputed and then the output from a filter rotated by any angle 9q can be found
by interpolating between them. When a large number of rotations of a template filter
is required, the above scheme can lead to substantial savings in both computational
cost (time) and memory requirements (space).
The necessary condition for a function j(i\9) to be steerable is that f(r,9) is
bandlimited in its polar angle:
N
f(r.9) = £ a n (r)e> n6
n=-N
it.17)
This can be verified by inserting (4.17) into (4.16) and by assuming, for con-
venience, that fi{r,9) = f(r,9 — 9,):
i
a n (r)e- jn6 ° =Y / c i (0o)a n {r)e- jnBi ,
(4.18)
i=i
where {#,-} is a set of user defined angles and n £ [— A", 0]. 4 Since only nonzero
coefficients a n (r) are of interest, both sides of (4.18) can be divided by a n (r). This
yields a matrix equation from which interpolating functions c 4 (# ) can be determined:
J6o
.3N60
=Jd
JNt
cJV2
o]N6 2
■ J" i
JNBj
ci(0o)
c 2 (0o)
ci(9 )
(4.19)
For coefficients a n = the rows corresponding to each n were removed from
the matrix formulation shown in (4.19). For this system to have a solution, the angles
{9;} must be chosen such that the columns of the matrix are linearly independent.
In [12] they proved that the minimum number of basis functions /,(;'. 9) needed
to steer /(r, 9) according to (4.16) is equal to the number of nonzero coefficients a n (r)
in the Fourier series expansion (4.17).
To conclude this brief description of steerability, let us only remark that func-
tions which are not steerable (i.e., do not have a finite number of terms in (4.17))
4 Note that the constraints are the same for n £ [-N, -1] and n£[l, .V], so that a subset of all
possible values for n € [-N, N] can be taken.
38
can be approximated with steerable functions (a singular value decomposition was
employed for approximating such functions efficiently [31]), and that the technique
of expressing transformed versions of a function as linear combinations of a fixed set
of basis functions is not limited to rotations (extensions to translations [39], scalings
[31, 39], and general transformations [14] have been reported).
4.3 Steerable Dyadic Wavelet Transform
Returning to the question from the beginning of Section 4.2, the answer seems
obvious: one needs to construct a steerable analog to the one-dimensional transform
from Chapter 3. Steerable transforms are nothing new — quite a few [13, 17, 38, 39]
have been devised, some of them [17, 38] exactly for the computation of directional
derivatives. Here, we are not interested in any directional derivatives: we want to
use derivatives of central B-splines which, as the order of B-splines increases, tend to
derivatives of a Gaussian.
We define a steerable dyadic wavelet transform of a function s(x, y) G L 2 {R 2 )
at a scale 2 m , m 6 Z, as [21]
Wi a s(x,y) = 8*tj>l(x,y), (4.20)
where ^(x,y) denotes ^ m (x,y) rotated by 9 U >:\„{x.y) = 2" 2m 0(2-"V,2- m y),
0(x,y) is a steerable wavelet that can be steered with I basis functions, and t = ^j-k
with i e {1,2,...,/}.
Analogously to the one-dimensional case (cf. Section 3.1) we require the two-
dimensional Fourier plane to be covered by the dyadic dilations of ^ i {2 m u x ,T n u y ):
there must exist ,4 3 > and B 3 < oo such that
4» < £ El^'(2 m ^,2 m u;,)| 2 < i?3 (4.21)
m=— oo i=l
is satisfied almost everywhere.
39
If (nonunique) reconstructing functions \^(:r,y) are chosen such that their
Fourier transforms satisfy
f) E ( / > '(2 m u; x ,2 m ^)x ! (2 m ^,2 m u; y ) = l, (4.22)
m=— oo t=l
the function s(x, y) may be reconstructed from its steerable dyadic wavelet transform
by
oo /
»{*,V)= E EwL**x\n(*,v)< (4-23)
m=— oo i=l
where \ l m {x,y) denotes \ m {x,y) rotated by 6, and \ m (x.y) — 2~ 2m \(2 _m .r, 2~ m y).
To derive an algorithm for the fast computation of the transform, we, similar
to (3.3), introduce two smoothing functions such that
00 I
i{L> x ,u> y )<p{u x ,u> y ) = £ ^^'(2 m ^,2 m a Jy )i'(2 m u; x ,2 m u; y ). (4.24)
We choose wavelets that are steerable analogs to the one-dimensional wavelets
from Section 3.1: 5
il>(u r ,ue) = (ju} r cos(u e )) d (— s^J , (4-25)
where u T = «/u;J + ^ and u^ = arg(cj :r ,u; y ). These wavelets can be steered with d+\
basis functions (i.e., / in (4.16) is equal to d+1).
Choosing
4>(2u r ) - H st (u r )<j)(u) r ), (4.26)
^ % {u} r ,u 6 ) = G st {uj r ,u: 9 - 9i) 6{u: r ). (4.27)
ip(2u T ) = L st (u> r )ip(u r ), (4.28)
and
X % (w r ,u e ) = K at (L> r ,u> g - Oi)<p{u) r ), (4.29)
This choice can be viewed as an extension of the transform presented in [20, 21, 24].
10
and using (4.26) through (4.29) with (4.24) computed for the finest two scales, we
obtain
^G a t(u; r ,u>4, - O^Kstiur,^ - 0i) + H 3t (ur)L s t{u T ) = 1. (4.30)
Setting <£(u> r ) = /? p (uv), and employing (3.5) and (4.25) with (4.26) and (4.27),
we find that
H«(u r ) = H-fa) (4.31)
and
G at {u r ,u e ) = (co S (u e )) d G- 3 {u> r ), (4.32)
where H(u>) and G(u) are specified by (3.9) and (3.15), respectively.
By inserting (4.31) and (4.32) into (4.30), the missing two filters can be chosen
as
L«(u r ) = L a (u r ) (4.33)
and
K at (u r ,U0) = — (cos{u e )) d K s {u r ), (4.34)
where L{^) and K(u)) are given by (3.16) and (3.17). respectively, and Cd =
Ef = i(cos(u;-^)) 2d .
4.4 Multiscale Spline Derivative-Based Transform
Let us pause here for a brief assessment of the two-dimensional steerable trans-
form derived so far. We have chosen steerable wavelets (4.25) which are equal to d-th
order derivatives of circularly symmetric spline functions in the direction of x-axis
(note that knots for these splines are circles) and we have laid a foundation for fil-
ter bank implementations in (4.30). An obvious shortcoming of this scheme is the
fact that none of the filter kernels obtained from (4.31) through ( 1.31) is compactly
supported on the rectangular grid. For implementations using digital filters, one is
therefore forced to approximate these frequency responses, and by doing so. (4.30)
II
may not hold anymore. Filters in filter bank implementations of steerable pyramids
described in [17, 38, 39], for example, were designed by using various techniques
for approximating desired frequency responses. None of the reported filter banks
achieved perfect reconstruction and all filter kernels were nonseparable. Here, we
will take a different approach. We will approximate the wavelets in (4.25) in a way
that will lead to x-y separable filters in a perfect reconstruction filter bank imple-
ment at ion of the transform such that the quality of approximation will improve by
increasing the order of B-splines.
Let us approximate wavelets from (4.25) with
, v d (3 p+ d(x)
nx,y) = — -jp — i3 P +d(y)- (4.35)
Based on the fact that B-splines tend to a Gaussian as their order increases,
it is easy to see that both wavelets (4.25) and (4.35) converge to the same functions
(i.e., d-ih order derivatives of the normalized Gaussian in the direction of .r-axis) as
p — * 'CO.
In order to steer wavelets tj){x,y) given by (4.35) (note that steering will
be only approximate, since these wavelets are not steerable), we need to find basis
functions that will approximately steer ij)(x,y). Computing rotations, as we did in
(4.20), is not practical here, because arbitrary rotations of (4.35) cannot be expressed
exactly in terms of central B-spline functions from Chapter 2. Instead, we take
advantage of the property of circularly symmetric functions that rotations of their
directional derivatives are equal to directional derivatives in rotated directions:
d d Q c (x,y)\ d d g c (x,y)
Tin
dn d dni
"0
where 1Z 6q stands for rotation by 6 , dec £' v) = nVg c (x,y), g c (x,y) is a circularly
symmetric function, and ng denotes vector n = (cos 0, sin 6) rotated by 6 .
Let us choose
g(x,y) = 8 p+d (x)/3 p+d {y),
12
which is approximately circularly symmetric function for higher order splines. A
rotation of %l>(x,y) = f^J y from (4.35) by 6q can therefore be approximated by
where n = (cos $o, sin #o) = (n x ,n y ).
Note that in case of Gaussian, which is both x-y separable and circularly
symmetric, (4.36) becomes exact (e.g., for g(x,y) = e~^ x + y ~>, o = —0, and d =
{2,4}, we obtain, up to a scaling factor, x-y separable basis set for the second and
fourth derivative of Gaussian from Tables III and VII of [12]).
Having found a set of basis functions (4.36) that approximately steer wavelets
(4.35), we want to construct a transform such that Equations (4.20) through (4.24)
will be valid (superscript i must be viewed now as an index, rather than rotation by
0i). In frequency domain, we can express basis functions from (4.36) as
4> i+1 (u x ,u; y ) = G±?(u x )GL,(Ljy)P P+i {u x )P P+i - i (u y ), £ = 0,1 d, (4.37)
where G d (u) is given by (3.15) and G°(w) = 1.
Choosing appropriate x % (u x ,u> y ) to obtain a relation needed for the filter bank
implementation of the transform is more complicated than in one dimension. Since we
could not find a general solution for arbitrary d, we solve each case separately. Below,
we present solutions for the first four orders. When d < 2, we impose ^(w^u,,) =
<f>(u> x ,u}y) = /3 p (u> x )fl p (u}y), a constraint analogous to the one from Section 3.1.
For d = 1 , we write similar to [28]
x\u>*,Uy) = K]{u x )T l {uj y )i3 p (uj x )fip-\{uy). (4.38)
\ 2 (u x ,^ y ) = r 1 (w x )j£(<«;J&_ 1 (w x )&(a; y ), (4.39)
where K d (u>) and Ti(u) are given by (3.17) and (4.11), respectively.
13
Computing (4.24) for the finest two scales and inserting (3.5). (4.37), (4.38),
and (4.39) yields a relation between frequency responses
G\u x )K\u> x )Ti{u>y) + TiK^Vv^Vv) + \H(u t )H{uy)\ 2 = 1.
For d = 2, we choose
\\u; r ,uj y ) = K*(u x )T 2 (uJy)Pp(u x )P P -2(u y ), (4.40)
X 2 {u x ,u v ) = K](u x )K}((jj y )P P - 1 (u g )P P -.i(u y ), (4.41)
X 3 (w*,w y ) = T 2 (u; ;r )A' s 2 (u;j / )^ p _2(u; ;c )^ p (LJ y ), (4.42)
w
here
r a (w) = |#(u,)| 2 . (4.43)
Using (3.5), (4.37). and (4.40) through (4.42) with (4.24) results in
G 2 (u x )K 2 (u x )T 2 (u y ) + G\u; x )K\u x )G l {w y )K\u y ) + T 2 (u x )G 2 (u y )K 2 (u y )+
+\H{u x )H(u> y )\ 2 = l.
For orders d > 2, we require <f>(u> x ,u y ) = J3 p (u} x )P P {u> y ) and ip(u x ,u> y ) = (f>(u} x )<p(uj y ),
where fi{u>) is specified by (3.7) and (3.16).
For d = 3. we choose reconstructing functions
^(u^Uy) = K*(w x )<p(u x )<p- 3 (L,j y ), (4.44)
X 2 (w„w y ) = -A' s 2 (^)A' s 1 (cj y )V3(u; r )V3(^)^-i(^ x )^-2(^), (4-45)
fK,w s ) = -A s 1 (^)A 2 (u; y )V3(a J;r )l / 3K)^-2(^)^-iK). (4.46)
X*(ws,w y ) = K*(u y )<p- 3 (u x )<f(u) y ), ( 1.47)
where
7 3 (w) = -^(1 - \H(lo)\ 2 ), (4.48)
and £_,-(a>) € L l (R) denotes a function such that <p(uj) = /3,_i(u;)<£_ t -(a;), i € iV.
II
Employing (4.37), (3.5), (3.7), and (4.44) through (4.47) with (4.24) yields a
relation
-G\u, x )K\u; x )V 3 (u; x )G\u: y )K 2 (u y )V 3 (u Jy ) + G 3 MK 3 M+
+H(u x )L(u> x )H(u y )L(u y ) = 1,
where L(*c) is specified by (3.16).
For d = 4, our choices are
X l (e>z,w y ) = Kt(oj x )T 2 {uj y )^(u) x )^- 4 M, (4.49)
X 2 (wx,u y ) = K 3 (u x )Kl(u y )<p-i(u x )<p- 3 (u> y ), (4.50)
X 3 (w*» w i») = - K Uvx)Kl(uy)V 4 (u x )V4(u y )<p-2(ux)<P-2(vv), (4.51)
X 4 (u x ,u y ) = A' s 1 (^ x )A' s 3 (u.' y )^_ 3 (u; J .)^_i(^ !/ ), (4.52)
X B (w t ,Uy) = T 2 (u> x )K 4 (u> y )<p- 4 (u x )ip(u v ), (4.53)
where
V 4 (u x ) = 1 - |tf(u;)| 2 . (4.54)
Using the above functions (4.49) through (4.53), (4.37), (3.5), and (3.7) in
(4.24) computed for the finest two scales gives
G 4 (u x )K 4 (u x )T 2 (u y ) + G 3 (u; a .)A' 3 (-.v)G 1 (u; J ,)A' 1 (u; J ,)-
-G 2 {u x )K 2 (u x )V 4 (u; x )G 2 {u y )K 2 {u; y )V 4 (u y ) + G\u^
+T 2 {uj x )G 4 {u; y )K 4 {uj y ) + H(u x )L{u x )H(u y )L(u> y ) = 1.
Here, we have even more freedom for choosing the reconstructing functions
than in one dimension. The above functions for d = {2,3,4} were found by trying
to imitate the one-dimensional transform from Chapter 3 as much as possible. All
decomposition filters G d (m) were first paired with corresponding reconstruction filters
45
K d (uj), and then all other potential digital filters were specified as polynomials in
H- a (u). We inserted thus specified filters into a relation between their frequency
responses and solved for the unknown polynomial coefficients. Since we allowed more
filters with higher-degree polynomials than expected in the solution, the resulting
system of linear equations was underdetermined. This allowed enough freedom for
removal of undesired digital filters and for balance between degrees of polynomials.
The described procedure for determination of reconstructing functions and
filters involves quite a lot of heuristics to obtain the appropriate solution from the
underdetermined linear system. Unfortunately, we are not aware of any systematic
way (aside from numerical optimization, which may be pretty cumbersome) to obtain
solutions comparable to the ones above.
Next, we will derive a filter bank implementation of the transform. As in
Section 4.1, we assume a bandlimited input signal: s(uj x ,u} y ) = for |o; z | > it or
|a> y | > 7r. Using Shannon's sampling theorem in two dimensions [16] with (4.20) and
basis functions from (4.37), we can write
/OO /"CO ^ -~
J2 J2 8 (**i »v)mhc(** - i x )ainc(t y - i y )-
-OO J — OO ■
l x = — OO ly = — OO
OO OO
Y. g d -7( m ?){3p+i(x -t x -m x ) Y. g- s {m y )f3 p+d -i{y -ty-m y )dt x dt y ,
m x = — oo m,u=—oo
where i = 0, 1, ...,(/ as in (4.37).
Again, we approximate sine functions with r-order cardinal splines, then use
(2.4), and get
x — n x ,y — fly
m—\
■B p+r+d . i+1 (u; y )G d - i (2 m u; x )GU2 m u y ) J] H p _ + a i (2 n u x )H^ d - i (2 n u; y ). (4.55)
71=0
Using (4.55) with an approximation 5 p+r+t+1 (w) ~ B p+r (uj)Bi(uj), we can ob-
tain a filter bank implementation of the transform decomposition. The reconstruction
part follows from (4.24), (4.37), and reconstructing functions for distinct values of d.
46
Figure 4.2 shows filter bank implementations of a multiscale spline derivative-based
transform for d — {1,2,3,4}. For d = 1, we recognize (except for the prefiltering
and postnltering) the filter bank implementation of a two-dimensional discrete dyadic
wavelet transform from [28]. For d = 2, however, our transform differs from the filter
bank presented in [22] (i.e., the corresponding transform described in Section 4.1):
second derivative is computed only in the directions of x and y-axis in [18, 22], which
is not enough for steering. Although not particularly appropriate for orientation anal-
ysis, such a scheme may still, as we have seen in Section 4.1, efficiently approximate
Laplacian of Gaussian across dyadic scales.
A transform similar to the one described in this section, was presented in
[17, 38, 39]. Their filter bank implementation is less redundant (downsampling is used
on the lowpass branch, while simultaneously keeping aliasing negligible by using a
filter with insignificant amount of energy for \uj r \ > |) and uses reconstruction filters
with same magnitude frequency responses as the decomposition ones — a possible
advantage for certain applications. They have, on the other hand, little control over
the function from which derivatives are computed (to obtain a d-th derivative, they
multiply a circularly symmetric filter by (-jcos6) d with all filters being obtained
by recursive minimization of a weighted combination of constraints), filter bank does
not have perfect reconstruction, and none of the filters is x-y separable.
IT
B> x ) B> y ) -JJ^hK) B p+r+1 (co y )
(a)
B p+r+1 (ro x ) B p+r+1 (co y ) B r (o) x ) B r (co y )
(b)
G 1 s (2 m co x )
G 1 s (2 m co v )
, sV - w y ,
H s (2 m (o x )H,(2V)
(c)
G 2 (2 m co x )
Qs(2 m a) x ) G 1 s (2 m co y )
G 2 (2 m co Y )
H,(2 m co x )H,(2V)
(e)
Kj2 m a> x ) T,(2 m co v )
"1 \ TT-l/^m
1,(2©,) K;(2 m co v )
L s (2 m o) x ) L s (2 m co y )
(d)
K 2 (2 m (o x )T 2 (2 m co v )
•'x/ ^V* - "'y/
K;(2 m co x )K;(2 m o) y )
T 2 (2 m co x ) K 2 (2 m co v )
L s (2 m co x ) L s (2 m co y )
(f)
Figure 4.2: Filter bank implementation of a multiscale spline derivative-based trans-
form for m G [0, M - 1]: (a) Prefiltering, (b) postfiltering, (c) decomposition and
(d) reconstruction modules for d = 1, and (e) decomposition and (f) reconstruction
modules for d = 2.
is
G 3 s (2 m co x ) B^V)
G 2 (2 m (o x )G.;(2 m (o y )
Gl(2 m co x ) G 2 (2 m co Y )
B 2 (2 m co x )Q 3 s (2 m (o y )-
H s (2 m (o x )H,(2 m co y )
(g)
K^coJ EfcV)
-K 2 (2 m (o x )V 3 (2 m co x )K;(2 m co,)V3(2 m (o y )
K;(2 m (o x )V 3 (2 m (o x )K 2 (2 ni co y )V3(2 m co y )
^(2 m co x )K 3 (2 m co y )
L s (2 m co x ) L S (2V)
(h)
G 4 (2 m (0 x ) B,(2 m co y )
G^X^coJB^V)
G 2 (2 m co x ) G 2 (2 m (o v )
G^coJ B 2 (2 m a) x ) C&V)
B 3 (2 m (o x ) G 4 (2 m co y )
H :S (2 m o) x )H s (23)
(i)
K 4 (2 m (o x )B 3 , (2 m (o y )T 2 (2V)
K 3 (2X)B,\2V)K:(2 m co y )
K 2 (2 m (o x )V 4 (2 m (o x )K 2 (2 m (o y )V 4 (2 m (o y )
R 2 , (2 m a) x )K;(2 m co x )K 3 (2 m (o y )
B 3 1 (2 m (o x )T 2 (2 m (o x )K 4 (2 m (0 y )
1 LsPX) L s (2 m co y )
0)
Figure 4.2: Continued: (g) Decomposition and (h) reconstruction modules for d = 3,
and (i) decomposition and (j) reconstruction modules for d = 4. Decomposition
modules are recursively applied at the locations of the filled circles.
CHAPTER 5
IMPLEMENTATION ISSUES
Except for the steerable dyadic wavelet transform, all transforms presented in
Chapters 3 and 4 can be implemented as filter banks consisting of one-dimensional fil-
ters only. In this chapter, we show how to take advantage of symmetry/antisymmetry
of filters when combined with a mirror extended input signal.
5.1 Finite Impulse Response Filters
Since all two-dimensional filters used in the filter bank implementations of
the transforms from Chapter 4 are x-y separable, we will begin this section with a
detailed description of FIR filter implementations for the one-dimensional discrete
dyadic wavelet transform implementation from Figure 3.1. The extension to two
dimensions will then be straightforward.
Let us refer to filters applied at scale 2 m as filters at level m+1, and let filters
at level 1 (Equations (3.9), (3.15) through (3.17), (4.11), (4.43), (4.48), and (4.54))
be called "original filters,'' to distinguish them from their upsampled versions. As
an input to the filter bank from Figure 3.1, we consider a real signal s(n) G 1 2 {Z).
n G [0, N - 1]. Depending on the length of each filter impulse response, filtering an
input signal may be computed either by multiplying the discrete Fourier transforms
of the two sequences or by circularly convolving s(n) with a filter's impulse response. 1
Of course, such a periodically extended signal may change abruptly at the boundaries
causing artifacts. A common remedy for such a problem is realized by constructing
x As is customary in image processing, we use circular, rather than linear, convolution.
1!)
50
a mirror extended signal [18]
s (n) = l S{ - n - l) if « € [-iV, -1]
mel j \ s(n) line [0,7V -1], ( ° j
where we chose the signal s me (n) to be supported in [— TV, N — 1]. It will become
evident shortly, that mirror extension is particularly elegant in conjunction with
symmetric/antisymmetric filters.
Let us first classify symmetric/antisymmetric real even-length signals into four
types [30]:
Type I f(n) = /(-»),
Type II /(n) = f(-n - 1),
Type III /(n) = -/(-«),
Type IV /( n ) = -/(-n-l),
where n € [-A r , A' - 1]. Note that for Type I signals the values at /(0) and f(—N)
are unique, and that for Type III signals the values at /(0) and f( — N) are equal to
zero.
Using properties of the Fourier transform it is easy to show that the convolu-
tion of symmetric/antisymmetric real signals results in a symmetric/antisymmetric
real signal. If a symmetric/antisymmetric real signal has an even length, then there
always exists an integer shift such that the shifted signal belongs to one of the above
types.
Now, we are ready to examine the filter bank implementation of a one-dimensional
discrete dyadic wavelet transform from Figure 3.1 with filters given by (3.9) and (3.15)
through (3.17) driven by a mirrored signal s me (n) at the input. Let the number of
levels M be restricted by
M < l + log 2 / ~* (5.2)
■'-'max -1-
51
where L max is the length of the longest original FIR filter impulse response.
Each FIR filter block in the filter bank consists of a filter and a circular shift
operator. Equation (5.2) guarantees that the length of the filter impulse response
does not exceed the length of the signal at any block.
Since our input signal s me (n) is of Type II and noninteger shifts at level 1 are
rounded to the nearest integer, it follows that a processed signal at any point in the
filter bank belongs to one of the types defined above. This means that filtering a
signal of length 2 A* can be reduced to filtering a signal of approximately one half of
its length. (For Types I and III, A + 1 samples are needed. However, for Type III
one needs to store only A — 1 values because zero values are always present at the
zeroth and (— JV)-th sample position).
Implementation is particularly simple for FIR filters designed with d even and
p odd. Filters are of Type I in this case, so the signal at any point of the filter bank
will be of Type II. An FIR filter block from the filter bank shown in Figure 3.1 can
therefore be implemented by
F s , m u(n) = f(0)u II (n) + Y,f(i)[u II (n-2 m i) + u JI (n + 2 m t)l n € [0, A r -1], (5.3)
where
u//(n) = <
u(-n-l) ifne[-f,-l]
u(n) ifn€[0,AT-l] (5.1)
u(2N-n-l) ifneJA,^-],
u(n) is an input signal to a block, /(??) is an impulse response of some original filter.
L is the length of the filter, and A* is the length of an input signal s(n) to the filter
bank. Implementation of filters b p (n) used for prefiltering and postfiltering represents
a special case of (5.3) with m = 0.
A filter bank with the above implementation of blocks and signal s(n) at the
input yields equivalent results as circular convolution for s me (n) as defined by (5.1).
In addition to requiring one half the amount of memory, the computational savings
52
over a circular convolution implementation of blocks are, depending on the original
filter length, three to four times fewer multiplications and one half as many additions.
A similar approach can be used for other filters. However, things get slightly
more complicated in this case, because the filters are not of the same type and the
signal components within the filter bank are of distinct types. As a consequence, an
implementation of blocks that use distinct original filters may not be the same, and
the implementation of blocks at level 1 may differ from the implementation of blocks
at other levels of analysis.
The decomposition blocks at level 1 can be implemented by
2 1
G- a , u(n) = J2 g{i)[u n {n - i - 1) - u n (n + i)], n G [1, JV - 1],
for d odd, (5.3) for d even,
H- s ,ou{n) = Y, h{i)[un{n-i-l) + u n (n + i)], n G [O.A],
t'=0
for p even, and (5.3) for p odd, where «//(/) is defined by (5.4), g(n) and h(n) are
impulse responses of the filters computed from (3.15) and (3.9), respectively, and L
is the length of the corresponding impulse response.
The output from a block G- 3 (u>) at level 1 is of Type III for d odd and of
Type II for d even, while the output from H_ a (u) at the same level is of Type I for
p even and of Type II for p odd.
The decomposition blocks at subsequent levels m G [l,Af — 1] can be imple-
mented by
G- a<m u(n) = £ g(i)[ui(n - 2 m (i + s)) - Ul (n + 2 m (i + s))}, n G [l,iV- 1],
;=0
for d odd and p even.
2 l
G- a , m u(n) = Y, 9(i)[un(n - 2 m {i + s)) - u n (n + 2 m (i + s))}, n G [0, N - 1],
f=0
53
L-\
2
for d and p odd.
F_ s>m «(n) = /(O)ur(n) + £ /COM" " 2™*') + "/(« + 2 "'0]- n € [0, AT], (5.5)
1=1
with f(n) = g(n) for d and p even,
H.,, m u(n) = J2 h(i)[ Ul (n - 2 m (i + a)) + u/(n + 2 TO (* + a))], n € [0, AT], (5.6)
i=0
for p even, and (5.3) for p odd, where
»(-??)
u/(n) == { u{n)
ifn€[-f,-l]
if n € [0, N\
(5.7
u(2N-n) \ine[N + l,m
The outputs from blocks G- S (2 m u:) are of Type III for d odd and p even, of
Type IV for d and p odd, and of Type I for d and p even, whereas the outputs from
H- S (2 m u)) are of Type I for p even and of Type II for p odd.
Next, the reconstruction blocks at level 1 can be implemented by
K s .ou{n) = J2 Hi)[uin(n -?' + !)- u ni (n + i)], n € [0, N - 1],
i=i
for d odd, (5.3) for r/ even,
L
2
L,, u(n) = £ '(»)[«/(» - »' + 1) + "/(» + *)]i n G [0, A r - 1],
! = 1
for p even, and (5.3) for p odd, where
uin(n) = <
-u(-n)
ifne[-f.
-11
if n =
u(n)
if n e [1,N
-1]
\(n = N
[5.8)
( -u(2tf-n) ifne [W+l,^],
«/(/?) is as defined by (5.7) and k(n) is an impulse response of the filter from (3.17).
Note that both outputs from blocks K a (u) and L a (u) are of Type II.
The reconstruction blocks at subsequent levels can be implemented by
i-l
2 *
K a<m u(n) = £ H* + !)[«///(» - 2 m (i + s)) - U//7 ( n + 2 m {i + .s))], n <= [0, AH,
!=0
54
for d odd and p even, (5.5) with f(n) = k{n) for d and p even,
A' s , m »(n) = £ fc(i + l)[u/v(n - 2 m (i + s))- u IV (n + 2 m (i + s))], n € [0. N - 1],
i=0
for d and p odd.
L sjn u(n) = #_ s , m u(n),
for p even, and (5.3) for p odd, where «///(/) is given by (5.8),
' -«(-„- 1) ifne[-f.-l]
u /v (n) = I n(n) if n € [O.N - 1]
, -«(2JV-n-l) ifnG[A',f-],
and H- S _ m u(n) is specified by (5.6). We observe that the outputs from blocks K s ('2 m uj)
and L s (2 m u;), m G [1. M — 1], are of Type I for p even, and of Type II for p odd.
When we compare the above implementation of blocks to circular convolution
driven by a mirrored signal s me (n) at the input, we observe that approximately
twofold less memory space, three to four times fewer multiplications and one half
as many additions are required. (For Type I signals an additional sample has to be
saved because two values are without a pair).
Until now. we have talked only about the one-dimensional case, whose filter
bank implementation is depicted in Figure 4.2. Two-dimensional transform filter
bank implementations (Figures 4.1 and 4.2) are not only comprised of x-y separable
filters solely, but also use all the filters from Section 3.1. Everything presented in
this section so far is therefore directly applicable to the two-dimensional case. Filters
which have nol been t reated yet (i.e.. /,(//). /_>(;/). r : ,,( n ). r,( n ). and filters b p (n ) from
the decomposition modules of Figure 4.2) can all be realized by (5.3) for p odd or
m = and (5.5) otherwise (/(??■) denotes an impulse response of any of the above
mentioned zero-phase filters in this case). 2
In case of filters i>3(n) and 1^(71), a slightly more efficient implementation can be obtained
by precomputing new filters k * i>3<n) and k * 1)4(71), and then implementing them by K s m u(n),
me [0,M- 1].
I 1 u- implementation presented in this section performs all operations in the
spatial domain, however, one could also implement the structures shown in Figures
3.1, 4.1, and 4.2 with an input signal s me (n) (Equation (5.1)) in the frequency domain.
For short filter impulse responses, such as those given in Tables 3.1. 3.2 and 3.3, the
spatial implementation described in this section is certainly more efficient. For long
filter impulse responses, however, filtering is faster if implemented in the frequency
domain. Additional details on alternative implementation strategies can be found in
[33].
5.2 Infinite Impulse Response Filters
Implementation of IIR filters b' 1 ^) which were introduced in Section 2.2 is
a bit more involved than the one encountered in the previous section. Fortunately,
the number of different cases is much smaller here: possible input to b^ln) in filter
banks from Figures 3.1, 4.1, and 4.2 is either of Type II or of Type I. 3 We will use
ideas and a few results from [43].
Let us first take a closer look at the system function B^iz). This function
can be written as a cascade of terms
E ^ - j^h^ - ,i- a --.) ( i- a -, - ^
which can be expressed in a parallel form as
-a / 1 1
-a 2 Vl -or 1 + Y^az
E ^ = ^-i 7— nr + i-^rr-0- <wo)
where a and - are poles of the causal and the anticausal filter, respectively.
The impulse response of this term can be written as
1 — cr
3 Note that symmetry types in this section slightly differ from those defined in Section 5.1: here,
mirror extended signals are periodically repeated, so that they stretch from — oo to oc.
.-,(,
We choose to implement E{z) in a cascade form and therefore extract the
difference equations from (5.9):
c + (n) = u(n) + ac + (n - 1) n = 1,2 V-l.
(5.11)
and
c(n) = a(c(n + l)-c + (n)) n = N-2,N-3,. . .,0, (5.12)
where u(n) denotes the input to the block, c + (n) is the output from the causal part,
and c(n) stands for the output from the block.
To solve (5.11) and (5.12) we need boundary conditions c + (0) and c(N — l).
Let us begin with filters fcr 1 (n) in filter bank implementations from Figures 3.1, 4.1,
4.2(a). 4.2(b). and Figures 4.2(h) and 4.2(j) with m = 0. We derive
■ ' v_1 a' +l + o
c+(0) = £ a-u„,(0 = «(0) + £ - T f^
:'=— oo i=0
and, using parallel form (5.10)
io
v27V
-</(/
u(0) + 5> <+1 ti(«), (5.13)
i'=0
c(N-r
—a
1 -a' 2
— a
1 -a 2
( c + (iV-i)+]T
-i n -V-, + Q .v+i +i
I
,2V
u(iY
i=0
N-l
[c + (N-l)+ Yl ^-^(i)),
i=.V — 1 — i'o
(5.14)
w
here
u 1Ip (n)
»//(") :
«//(;? mod (27V)) if n>0
»//(-(n + l) mod (2JV)) if n < 0,
u(n) if n e [0, A r - 1]
u(2JV-n- 1) if/? € [.V.2.V- 1],
iV is the length of an input signal to the filter bank, and ?o < N — l is selected such
that a la falls below a predefined precision threshold.
For IIR filters from Figures 4.2(h) and 4.2(j) with m G [1. M - 1] and p odd,
we get
o
c+ (°) = £ [a-**u IIp (i
t=— oo
oo N-l
Y.Y.
n=0 i=0
OO iV— 1 .. r -, N
.„. W--V ^^ I r »+l+2Nn 1 2.Y(n+l)-i |
= "(0) + EE [o^ 5 ^] + a-^r— u (i)
(5.15)
57
and
N-\
^(A f -i) = T^ T (^(-v-i) + EL{h
2?n
+
n—0 i'=0
Q
.Y+l + i + 2.Yn
}«(0, (5-16)
where
[«1u -
n
x if .r € Z
otherwise.
If N is a power of two and 2 m_1 < JV ((5.2) guarantees that the latter condition
is always true) (5.15) becomes
c+(o) = «(o)+ Y.
i+1
a 2 m
+
2-V-i
Q 2"
;=o
W
hile (5.16) can be written as
c(X-l) = -^(c + (N-l)+Y
— o ^~^
1 — Q 2™
N-i
Q 2 m
a «(0,
+
■ v tH'
Q 2 m
1 - Q2 m
-»(/)•
Finally, the boundary conditions for filters ^(rc) from Figures 4.2(h) and
4.2(j) with m E [1,M — 1] and p even are
c+(0)= £ [a-^u/pCi)
i=— oo
V^ f r 2ATn 1
^ j [a -5™-] i/(0)+ a^ 7-
+
n=0
N-l
+ 2JVn
2A r |7i + l)-
Q 2">
u(N)+
u(i)\
(5.17)
and
where
—a
c(7V) = - -(2c + (AT)-«(JV)),
1 -a
uj(n) =
«/p(«) = »/(H mod (2A*)).
' u(n) if n € [0. A"]
u(2JV-n) if n e [N + 1,2N-1],
and c~(N) = c + {N) with c~(n) denoting response of the anticausal filter from (5.10)
was used.
oS
Again, if N is a power of two and 2 m l < N, (5.17) can be simplified
+ m\ —
c + (0)
u(0) + ai^u(N) + T.l^{
Q2 r
+
2iV-i
}tl(l
1 2i£
1 — a 2™
Series in expressions for c + (0), c(N — 1). and c(N) for filters from Figures
4.2(h) and 4.2(j) with m G [1. M — 1] can be, similar to (5.13) and (5.14), truncated
according to the desired precision.
For orders p greater than three, we implement B~ 1 {z) as a cascade of terms
E{z) with different q's. Note that the output from block E(z) is always of the same
type as the input to it.
CHAPTER 6
APPLICATIONS
In this chapter, we present a. comparison between multiscale spline derivative-
based transform and traditional wavelet techniques on two sample applications.
6.1 Image Fusion
Image fusion combines particular aspects of information from the same imag-
ing modality or from distinct imaging modalities and can be used to improve the
reliability of a particular computational vision task or to provide a human observer
with a deeper insight about the nature of observed data. Whether it is combining
different sensors or extending the dynamic range of a single sensor, the goal is to
achieve more accurate inferences that can be achieved by a single sensor or a single
sensor setting.
The simplest method of fusing images is accomplished by computing their
average. Such a technique does combine features from input images in the fused
image, however, the contrast of the original features can be significantly reduced.
Among more sophisticated methods, multiscale and multiresolution analyses have
become particularly popular. Different pyramids [5, 41] and wavelet-based techniques
[7, 20, 26. 32] have been applied to this problem. Fusion is performed in the transform
space by computing local statistics across the decomposition scales. Typical size of
neighborhood is between a single pixel and 5 x 5 area with some loss of contrast
reported for the latter [5]. Some authors argue about of the particular criteria for
fusion [5, 26], while others propose a set of different combination rules to perform
image fusion in a variety of ways [7].
:»!)
60
Here, we are not trying to find the best criteria for fusion. We simply use
a "maximum magnitude" rule (i.e., at each position and scale of the transform the
corresponding transform coefficients are compared and those with greater magnitude
included for reconstruction) and compare the results obtained by traditional wavelet
techniques with the ones produced by the multiscale spline derivative- based trans-
form.
Figure 6. 1(a), (b) shows a pair of images with distinct (lower in Figure 6.1(a)
and upper in Figure 6.1(b)) areas blurred. Figure 6.1(c) demonstrates the ideal
result of image fusion obtained by manually cutting and pasting the upper part of
Figure 6.1(a) and the lower part of Figure 6.1(b). Image fusion resulting from the use
of discrete wavelet transforms is shown in Figure 6.1(d) for the orthogonal wavelet
DAUB12 [7] and in Figure 6.1(e) for the linear phase biorthogonal wavelet Bior6.8.
Figure 6.1(f) illustrates the result of image fusion with a multiscale spline derivative-
based transform being employed. Note that the image fused using the multiscale
spline derivative-based transform is virtually indistinguishable from the one obtained
manually. Both orthogonal and biorthogonal wavelet transform exhibit artifacts. On
the average, biorthogonal transform performs better than orthogonal due to the linear
phase of wavelets and filters, but still causes serious artifacts.
Basing the decision rule on an area rather than on a single pixel certainly
reduces artifacts of traditional wavelet methods, however, it reduces sharpness as
well. As noted by Burt and Kolczynski [5], "there is a tradeoff between sharpness
and shift invariance." Here, we argue the point that the proper choice of transform
can eliminate artifacts due to translation and rotation invariance altogether.
Another possible source of artifacts is misregistration. Low-level fusion al-
gorithms typically require input images that have already been properly aligned.
Misregistration causes artifacts [7], and it can be recognized from our discussion
1,1
in Section 1.1 that they are more obvious for translation and rotation-noninvariant
wavelet techniques than for the multiscale spline derivative-based transform.
6.2 Image Enhancement
Similar to the previous section, we will point out problems associated with
traditional wavelet analysis on an example. Below, we briefly present the problem of
contrast enhancement in digital mammography.
In mammography, early detection of breast cancer relies upon the ability to
distinguish between malignant and benign mammographic features. Contrast en-
hancement can make more obvious unseen or barely seen features of a mammogram
without requiring additional radiation. The complexity of mammographic images
and the subtlety of malignancies can present a challenge even to expert radiologists.
In addition to dealing with barely visible mammographic features, human observers
sometimes simply overlook abnormalities. Such mistakes affect the number of false-
negative cases considerably. Computer processing of digital mammograms can assist
radiologists to reach a correct diagnosis more consistently. A variety of computer
based techniques have been reported in almost three decades of research. The advent
of direct digital mammography devices has made digital image processing techniques
more attractive for screening.
Contrast enhancement in a wavelet transform framework can be achieved by
applying a (typically nonlinear) function to wavelet coefficients and then reconstruct-
ing an enhanced image with modified coefficients. We are not trying to find the best
method for a specific type of images or noise here, rather we want to illustrate pos-
sible sources of artifacts when orthogonal and biorthogonal wavelet transforms are
used.
Figure 6.2(a) shows an original mammographic image. Figure 6.2(b) depicts
the result of enhancement using piecewise linear enhancement function [11, 23] to
62
(c)
(d)
Figure 6.1: (a) An image with lower part blurred, (b) An image with upper part
blurred, (c) Fused image obtained by combining images from (a) and (b) man-
ually. Results of fusion performed by (d) discrete orthogonal wavelet transform
(DAUB12), (e) discrete biorthogonal wavelet transform (Bior6.8), and (f) multiscale
spline derivative-based transform.
63
modify wavelet transform coefficients, and Figure 6.2(c) shows the result of enhance-
ment using the multiscale spline derivative-based transform and the same enhance-
ment function as Figure 6.2(b). In this figure, we. more than trying to point out
better enhancement results in Figure 6.2(c). observe that Figure 6.2(b) shows fea-
tures (artifacts) that are not present in the original image (they are particularly
obvious in the area right from the mass in Figure 6.2(b)). Here, the goal of enhance-
ment is making visible features that are present in the original image but obstructed
by different structures within the same image. Generation of artificial features, while
not being rare when nonlinearly processing coefficients of the orthogonal or biorthog-
onal wavelet transform, is certainly not desirable for this particular application. For
more details on different enhancement strategies in mammography, please refer to
[11. 22, 23, 25].
Another problem related to enhancement and other kinds of processing of dig-
ital mammograms (and not unique to mammography) is that the same mammogram
could have been easily translated and/or rotated. From obvious reasons, it is not
desirable that the result of processing can be significantly affected by positioning
only — one more reason for a translation and rotation-invariant transform.
Figure 6.3 shows outputs from steerable basis filters applied to the image from
Figure 6.2(a) at scales {1.2.4}. By linearly combining these outputs analysis along
arbitrary direction can be performed (cf. Chapter 4). This enables, as illustrated in
Chapter 1. a rotation-invariant processing [21].
64
Figure 6.2: (a) An original mammographic image containing a spicular mass. An
enhanced image using (b) discrete biorthogonal wavelet transform (Bior6.8) and (c)
multiscale spline derivative-based transform.
(>.->
Figure 6.3: The magnitude of filter coefficients for G%'{u} T ,u<jy) at fc = {0, \t, \tc]
(three levels of analysis shown).
CHAPTER 7
CONCLUSION
In this thesis we constructed a new transform that does not introduce artifacts
due to translation and rotation invariance, which are inherent to traditional wavelet
analysis.
We extended the one-dimensional discrete dyadic wavelet transform to higher-
order derivatives and even-order spline functions and developed an improved initial-
ization procedure. Comparison to the originally employed initialization [28] showed
significantly better performance of our procedure for finer scales of analysis.
We developed several two-dimensional transforms. All of them were derived
with the goal of eliminating artifacts due to lack of translation and rotation in-
variance. We presented a two-dimensional discrete dyadic wavelet transform with a
first-derivative wavelet as an extension of the one originally proposed in [28], a two-
dimensional discrete dyadic wavelet transform with a second-derivative wavelet that
can approximate Laplacian of a. Gaussian, and a steerable dyadic wavelet transform
implemented in a near-perfect reconstruction filter bank which may be preferred when
there is a need for orientation processing along equally spaced angles. We derived
a multiscale spline derivative-based transform which uses x-y separable filters in a
perfect reconstruction filter bank and enables fast translation and rotation-invariant
directional analysis of images.
We empirically showed the presence of artifacts in image fusion and enhance-
ment applications when traditional wavelet methods were used. Here developed trans-
forms did not exhibit similar artifacts, and, in case of fusion, yielded nice results
without the use of postprocessing.
66
67
Future work will concentrate on a more thorough exploitation of the transform
space for the existing applications, and potential new applications (e.g., denoising). In
particular, we will try to use the causality property of a Gaussian kernel mentioned
in Chapter I to trace features across scales and process them according to their
multiscale behavior.
REFERENCES
1] A. Aldroubi, M. Unser, and M. Eden, "Cardinal spline filters: Stability and
convergence to the ideal sine interpolator," Signal Process., vol. 28, no. 2, pp.
127-138, 1992.
2] J. Babaud. A. P. Witkin, M. Baudin, and R. 0. Duda, "Uniqueness of the
Gaussian kernel for scale-space filtering," IEEE Trans. Pattern Anal. Mach.
Intel!., vol. 8. no. 1, pp. 26-33, 1986.
3] C. de Boor, A Practical Guide to Splines, Springer- Verlag, New York, NY, 1978.
4] P. J. Burt and E. H. Adelson, "The Laplacian pyramid as a compact image
code," IEEE Trans. Commun., vol. 31, no. 4, pp. 532-540, 1983.
5] P. J. Burt and R. J. Kolczynski, "Enhanced image capture through fusion," in
Proc. Int. Conf. Comput. Vision, Berlin, Germany, 1993, pp. 173-182.
6] J. Canny, "A computational approach to edge detection," IEEE Trans. Pattern
Anal. Mach. IntelL, vol. 8. no. 6, pp. 679-698, 1986.
7] L. J. Chipman, T. M. Orr, and L. N. Graham, "Wavelets and image fusion,"
in Proc. IEEE Int. Conf. Image Process., Washington, D.C., 1995, vol. 3, pp.
248-251.
8] C. K. Chui, An Introduction to Wavelets, Academic Press, San Diego, CA, 1992.
9] C. K. Chui, Ed., Wavelets: A Tutorial in Theory and Applications, Academic
Press. San Diego, CA, 1992.
[10] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, PA, 1992.
[\Y\ J. Fan and A. Laine, "Multiscale contrast enhancement and denoising in digital
radiographs," in Wavelets in Medicine and Biology, A. Aldroubi and M. Unser,
Eds.. CRC Press, Boca Raton, FL, 1996, pp. 163-189.
[12] W. T. Freeman and E. H. Adelson, "The design and use of steerable filters,"
IEEE Trans. Pattern Anal. Mach. IntelL, vol. 13, no. 9, pp. 891-906, 1991.
[13] H. Greenspan, S. Belongie, R. Goodman, P. Perona, S. Rakshit, and C. H.
Anderson, "Overcomplete steerable pyramid filters and rotation invariance," in
Proc. IEEE Comput. Soc. Conf. Comput. Vision Pattern Recognit., Seattle, WA,
1994, pp. 222-228.
[14] Y. Hel-Or and P. C. Teo, "Canonical decomposition of steerable functions," in
Proc. IEEE Comput. Soc. Conf. Comput. Vision Pattern Recognit., San Fran-
cisco, CA, 1996, pp. 809-816.
69
[15] M. Holschneider, R. Kronland-Martinet, J. Morlet, and Ph. Tchamitchian, "A
real-time algorithm for signal analysis with the help of the wavelet transform," in
Wavelets, Time-Frequency Methods and Phase Space, J. M. Combes, A. Gross-
mann, and Ph. Tchamitchian, Eds., Springer- Verlag, Berlin, Germany, 1989, pp.
286-297.
[16] A. J. Jerri, "The Shannon sampling theorem — its various extensions and appli-
cations: A tutorial review," Proc. IEEE, vol. 65, no. 11, pp. 1565-1596, 1977.
[17] A. Karasaridis and E. Simoncelli, "A filter design technique for steerable pyramid
image transforms," in Proc. IEEE Int. Conf Acoust. Speech Signal Process.,
Atlanta, GA, 1996, vol. 4, pp. 2389-2392.
[18] I. Koren and A. Laine, "A discrete dyadic wavelet transform for multidimensional
feature analysis," in Time-Frequency and Wavelet Transforms in Biomedical
Engineering, M. Akay, Ed.. IEEE Press, New York, NY, 1997.
[19] I. Koren, A. F. Laine, J. Fan, and F. J. Taylor, "Edge detection in echocardio-
graphic image sequences by 3-d multiscale analysis," in Proc. IEEE Int. Conf.
Mage Process., Austin, TX, 1994, vol. 1, pp. 288-292.
[20] I. Koren. A. Laine. and F. Taylor, "Image fusion using steerable dyadic wavelet
transform," in Proc. IEEE Int. Conf. Image Process., Washington, D.C., 1995,
vol. 3, pp. 232-235.
21}) I. Koren, A. Laine. F. Taylor, and M. Lewis, "Interactive wavelet processing and
techniques applied to digital mammography," in Proc. IEEE Int. Conf. Acoust.
Speech Signal Process., Atlanta, GA, 1996, vol, 3, pp. 1415-1418.
] A. Laine, J. Fan, and S. Schuler, "A framework for contrast enhancement by
dyadic wavelet analysis." in Digital Mammography, A. G. Gale, S. M. Astley,
D. R. Dance, and A. Y. Cairns, Eds., Elsevier, Amsterdam, The Netherlands,
1994, pp. 91-100.
A. Laine, J. Fan, and W. Yang, "Wavelets for contrast enhancement of digital
mammography," IEEE Eng. Med. Biol. Mag., vol. 14, no. 5, pp. 536-550, 1995.
[24] A. Laine, I. Koren, W. Yang, and F. Taylor, "A steerable dyadic wavelet
transform and interval wavelets for enhancement of digital mammography," in
Wavelet Applications II. Proc. SPIE, H. H. Szu, Ed., Orlando, FL, 1995, vol.
2/91, pp. 736-749.
[25] A. F. Laine, S. Schuler, J. Fan, and W. Huda, "Mammographic feature en-
hancement bv multiscale analysis," IEEE Trans. Med. Imaging, vol. 13. no. 4,
pp. 725-740, 1994.
[26] H. Li, B. S. Manjunath, and S. K. Mitra, "Multi-sensor image fusion using the
wavelet transform," in Proc. IEEE Int. Conf. Image Process., Austin, TX, 1994,
vol. 1, pp. 51-55.
[27] S. Mallat and W. L. Hwang, "Singularity detection and processing with
wavelets," IEEE Trans. Inf. Theory, vol. 38, no. 2, pp. 617-643, 1992.
[28] S. Mallat and S. Zhong, "Characterization of signals from multiscale edges,"
IEEE Trans. Pattern Anal. Mach. Inteli, vol. 14, no. 7, pp. 710-732, 1992.
71)
[29] D. Marr and E. Hildreth, "Theory of edge detection," Proc. R. Soc. London
Ser. B, vol. 207, no. 1167. pp. 187-217, 1980.
[30] A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, Prentice-
Hall, Englewood Cliffs, NJ, 1989.
[31] P. Perona. "Deformable kernels for early vision," IEEE Trans. Pattern Anal.
Mach. Intel!., vol. 17, no. 5, pp. 488-499, 1995.
[32] T. Ranchin, L. Wald, and M. Mangolini, "Efficient data fusion using wavelet
transform: the case of SPOT satellite images," in Matin nmtical Imaging:
Wavelet Applications in Signal and Image Processing. Proc. SPIE. A. F. Laine,
Ed., San Diego, CA, 1993, vol. 2034, pp. 171-178.
[33] 0. Rioul and P. Duhamel, "Fast algorithms for discrete and continuous wavelet
transforms," IEEE Trans. Inf. Theory, vol. 38, no. 2, pp. 569-586, 1992.
[34] I. J. Schoenberg, "Contributions to the problem of approximation of equidistant
data by analytic functions," Q. Appl. Math., vol. 4. no. 1. pp. 45-99, 112-141,
1946. '
[35] I. J. Schoenberg. "Cardinal interpolation and spline functions." ./. Approx.
Theory, vol. 2, no. 2, pp. 167-206, 1969.
[36] I. J. Schoenberg, "Notes on spline functions III: On the convergence of the
interpolating cardinal splines as their degree tends to infinity," Isr. J. Math..
vol. 16. no. 1, pp. 87-93, 1973.
[37] C. E. Shannon, "Communication in the presence of noise," Proc. IRE. vol. 37,
no. 1, pp. 10-21, 1949.
[38] E. P. Simoncelli and W. T. Freeman, "The steerable pyramid: A flexible archi-
tecture for multi-scale derivative computation," in Proc. IEEE Int. Conf. Image
Process., Washington, D.C., 1995, vol. 3, pp. 444-447.
[39] E. P. Simoncelli, W. T. Freeman, E. H. Adelson, and D. J. Heeger, "Shiftable
multiscale transforms," IEEE Trans. Inf. Theory, vol. 38, no. 2, pp. 587-607,
1992.
[40] M. J. T. Smith and T. P. Barnwell, III, "Exact reconstruction techniques for
tree-structured subband coders," IEEE Trans. Acoust. Speech Signal Process.,
vol. 34, no. 3, pp. 434-441, 1986.
[41] A. Toet, "Image fusion by a ratio of low-pass pyramid." Pattern Recognit. Lett.,
vol. 9. no. 4, pp. 245-253, 1989.
[42] M. Unser and A. Aldroubi, "A general sampling theory for nonideal acquisition
devices," IEEE Trans. Signal Process., vol. 42, no. 11, pp. 2915-2925, 1994.
[43] M. Unser, A. Aldroubi, and M. Eden, "Fast B-spline transforms for continuous
image representation and interpolation," IEEE Trans. Pattern Anal. Mach.
Intel!., vol. 13. no. 3, pp. 277-285, 1991.
[44] M. Unser, A. Aldroubi. and M. Eden, "On the asymptotic convergence of B-
spline wavelets to Gabor functions," IEEE Trans. Inf. Theory, vol. 38, no. 2,
pp. 864-872, 1992.
71
[45] M. Unser, A. Aldroubi, and M. Eden, "Polynomial spline signal approximations:
Filter design and asymptotic equivalence with Shannon's sampling theorem,"
IEEE Trans. Inf. Theory, vol. 38, no. 1, pp. 95-103, 1992.
[46] M. Unser, A. Aldroubi, and S. J. Schiff, "Fast implementation of the continuous
wavelet transform with integer scales," IEEE Trans. Signal Process., vol. 42, no.
12, pp. 3519-3523, 1994.
[47] M. Vrhel, C. Lee, and M. Unser, "Fast computation of the continuous wavelet
transform through oblique projections," in Proc. IEEE Int. Conf. Acoust. Speech
Signal Process., Atlanta, GA, 1996, vol. 3, pp. 1459-1462.
[48] A. Witkin. "Scale space filtering," in Proc. Int. Joint Conf. Artif. IntelL,
Karlsruhe, Germany, 1983, pp. 1019-1022.
BIOGRAPHICAL SKETCH
Iztok Koren earned the B.S. and M.S. degrees in electrical engineering from the
University of Ljubljana, Slovenia, in 1987 and 1991, respectively. He will receive the
Ph.D. degree in electrical and computer engineering from the University of Florida,
Gainesville, in 1996.
His research interests include image processing, digital signal processing, and
wavelets.
He is a member of Eta Kappa Nu, Tau Beta Pi, and the IEEE.
7-2
I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate, in scope and
quality, as a dissertation for the degree of Doctor of Philosophy.
h/4/ /f4f/fr^
,Fred J/Taylor/, Chairman
Profess'or of/Electrical and Computer
Engineering
I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate, in scope and
quality, as a dissertation for the degree of Doctor of Philosophy.
6^/iU. yj , -<^£Ui>
Andrew F. Laine , Cochairman
Associate Professor of Computer and
Information Science and Engineering
I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and-^-fylly adequate, in scope and
quality, as a dissertation for the degree of Doctor
/
Jose C/
olr of Ele^rical and Computer
ngmeering
I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate, in scope and
quality, as a dissertation for the degree of DocJ*^ of Philosophy.
'John M. M. Anderson
Assistant Professor of Electrical and
Computer Engineering
I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate, in scope and
quality, as a dissertation for the degree of Doyetor of Philosophy.
■pai
6/lA
?mit Sigmon
Associate Professor of Mathematics
This dissertation was submitted to the Graduate Faculty of the College of
Engineering and to the Graduate School and was accepted as partial fulfillment of
the requirements for the degree of Doctor of Philosophy.
December 1996
V^ Winfred M. Phillips
Dean, College of Engineering
Karen A. Holbrook
Dean, Graduate School
L
LD
1780
19S£
.Km
UNIVERSITY OF FLORIDA
3 1262 08555 0696
TVNews
OpenLibrary